Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for reading and writing restart files.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_io
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE basis_set_types, ONLY: gto_basis_set_type
16 : USE cell_types, ONLY: cell_type
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_csr_create_from_dbcsr, &
19 : dbcsr_csr_dbcsr_blkrow_dist, dbcsr_csr_destroy, dbcsr_csr_type, dbcsr_csr_write, &
20 : dbcsr_desymmetrize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_has_symmetry, dbcsr_release, &
21 : dbcsr_type
22 : USE cp_files, ONLY: close_file,&
23 : open_file
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_get_default_io_unit,&
26 : cp_logger_type
27 : USE cp_output_handling, ONLY: cp_p_file,&
28 : cp_print_key_finished_output,&
29 : cp_print_key_should_output,&
30 : cp_print_key_unit_nr
31 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
32 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
33 : section_vals_type,&
34 : section_vals_val_get
35 : USE kinds, ONLY: default_path_length,&
36 : default_string_length,&
37 : dp
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE pao_input, ONLY: id2str
40 : USE pao_param, ONLY: pao_param_count
41 : USE pao_types, ONLY: pao_env_type
42 : USE particle_types, ONLY: particle_type
43 : USE physcon, ONLY: angstrom
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_kind_types, ONLY: get_qs_kind,&
47 : pao_potential_type,&
48 : qs_kind_type
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_io'
56 :
57 : PUBLIC :: pao_read_restart, pao_write_restart
58 : PUBLIC :: pao_read_raw, pao_kinds_ensure_equal
59 : PUBLIC :: pao_ioblock_type, pao_iokind_type
60 : PUBLIC :: pao_write_ks_matrix_csr, pao_write_s_matrix_csr
61 : PUBLIC :: pao_write_hcore_matrix_csr, pao_write_p_matrix_csr
62 :
63 : ! data types used by pao_read_raw()
64 : TYPE pao_ioblock_type
65 : REAL(dp), DIMENSION(:, :), ALLOCATABLE :: p
66 : END TYPE pao_ioblock_type
67 :
68 : TYPE pao_iokind_type
69 : CHARACTER(LEN=default_string_length) :: name = ""
70 : INTEGER :: z = -1
71 : CHARACTER(LEN=default_string_length) :: prim_basis_name = ""
72 : INTEGER :: prim_basis_size = -1
73 : INTEGER :: pao_basis_size = -1
74 : INTEGER :: nparams = -1
75 : TYPE(pao_potential_type), ALLOCATABLE, DIMENSION(:) :: pao_potentials
76 : END TYPE pao_iokind_type
77 :
78 : INTEGER, PARAMETER, PRIVATE :: file_format_version = 4
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief Reads restart file
84 : !> \param pao ...
85 : !> \param qs_env ...
86 : ! **************************************************************************************************
87 8 : SUBROUTINE pao_read_restart(pao, qs_env)
88 : TYPE(pao_env_type), POINTER :: pao
89 : TYPE(qs_environment_type), POINTER :: qs_env
90 :
91 : CHARACTER(LEN=default_string_length) :: param
92 : INTEGER :: iatom, ikind, natoms
93 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2kind
94 8 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
95 : LOGICAL :: found
96 : REAL(dp) :: diff
97 8 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hmat, positions
98 8 : REAL(dp), DIMENSION(:, :), POINTER :: block_X, buffer
99 : TYPE(cell_type), POINTER :: cell
100 : TYPE(mp_para_env_type), POINTER :: para_env
101 8 : TYPE(pao_ioblock_type), ALLOCATABLE, DIMENSION(:) :: xblocks
102 8 : TYPE(pao_iokind_type), ALLOCATABLE, DIMENSION(:) :: kinds
103 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
104 :
105 0 : CPASSERT(LEN_TRIM(pao%restart_file) > 0)
106 8 : IF (pao%iw > 0) WRITE (pao%iw, '(A,A)') " PAO| Reading matrix_X from restart file: ", TRIM(pao%restart_file)
107 :
108 : CALL get_qs_env(qs_env, &
109 : para_env=para_env, &
110 : natom=natoms, &
111 : cell=cell, &
112 8 : particle_set=particle_set)
113 :
114 : ! read and check restart file on first rank only
115 8 : IF (para_env%is_source()) THEN
116 4 : CALL pao_read_raw(pao%restart_file, param, hmat, kinds, atom2kind, positions, xblocks)
117 :
118 : ! check cell
119 52 : IF (MAXVAL(ABS(hmat - cell%hmat)) > 1e-10) THEN
120 0 : CPWARN("Restarting from different cell")
121 : END IF
122 :
123 : ! check parametrization
124 4 : IF (TRIM(param) /= TRIM(ADJUSTL(id2str(pao%parameterization)))) THEN
125 4 : CPABORT("Restart PAO parametrization does not match")
126 : END IF
127 :
128 : ! check kinds
129 11 : DO ikind = 1, SIZE(kinds)
130 11 : CALL pao_kinds_ensure_equal(pao, qs_env, ikind, kinds(ikind))
131 : END DO
132 :
133 : ! check number of atoms
134 4 : IF (SIZE(positions, 1) /= natoms) THEN
135 0 : CPABORT("Number of atoms do not match")
136 : END IF
137 :
138 : ! check atom2kind
139 15 : DO iatom = 1, natoms
140 15 : IF (atom2kind(iatom) /= particle_set(iatom)%atomic_kind%kind_number) THEN
141 0 : CPABORT("Restart atomic kinds do not match.")
142 : END IF
143 : END DO
144 :
145 : ! check positions, warning only
146 4 : diff = 0.0_dp
147 15 : DO iatom = 1, natoms
148 48 : diff = MAX(diff, MAXVAL(ABS(positions(iatom, :) - particle_set(iatom)%r)))
149 : END DO
150 4 : CPWARN_IF(diff > 1e-10, "Restarting from different atom positions")
151 :
152 : END IF
153 :
154 : ! scatter xblocks across ranks to fill pao%matrix_X
155 : ! this could probably be done more efficiently
156 8 : CALL dbcsr_get_info(pao%matrix_X, row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
157 30 : DO iatom = 1, natoms
158 88 : ALLOCATE (buffer(row_blk_sizes(iatom), col_blk_sizes(iatom)))
159 22 : IF (para_env%is_source()) THEN
160 11 : CPASSERT(row_blk_sizes(iatom) == SIZE(xblocks(iatom)%p, 1))
161 11 : CPASSERT(col_blk_sizes(iatom) == SIZE(xblocks(iatom)%p, 2))
162 193 : buffer = xblocks(iatom)%p
163 : END IF
164 750 : CALL para_env%bcast(buffer)
165 22 : CALL dbcsr_get_block_p(matrix=pao%matrix_X, row=iatom, col=iatom, block=block_X, found=found)
166 22 : IF (ASSOCIATED(block_X)) THEN
167 375 : block_X = buffer
168 : END IF
169 52 : DEALLOCATE (buffer)
170 : END DO
171 :
172 : ! ALLOCATABLEs deallocate themselves
173 :
174 34 : END SUBROUTINE pao_read_restart
175 :
176 : ! **************************************************************************************************
177 : !> \brief Reads a restart file into temporary datastructures
178 : !> \param filename ...
179 : !> \param param ...
180 : !> \param hmat ...
181 : !> \param kinds ...
182 : !> \param atom2kind ...
183 : !> \param positions ...
184 : !> \param xblocks ...
185 : !> \param ml_range ...
186 : ! **************************************************************************************************
187 21 : SUBROUTINE pao_read_raw(filename, param, hmat, kinds, atom2kind, positions, xblocks, ml_range)
188 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
189 : CHARACTER(LEN=default_string_length), INTENT(OUT) :: param
190 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hmat
191 : TYPE(pao_iokind_type), ALLOCATABLE, DIMENSION(:) :: kinds
192 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2kind
193 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: positions
194 : TYPE(pao_ioblock_type), ALLOCATABLE, DIMENSION(:) :: xblocks
195 : INTEGER, DIMENSION(2), INTENT(OUT), OPTIONAL :: ml_range
196 :
197 : CHARACTER(LEN=default_string_length) :: label, str_in
198 : INTEGER :: i1, i2, iatom, ikind, ipot, natoms, &
199 : nkinds, nparams, unit_nr, xblocks_read
200 : REAL(dp) :: r1, r2
201 : REAL(dp), DIMENSION(3) :: pos_in
202 : REAL(dp), DIMENSION(3, 3) :: hmat_angstrom
203 :
204 21 : CPASSERT(.NOT. ALLOCATED(hmat))
205 21 : CPASSERT(.NOT. ALLOCATED(kinds))
206 21 : CPASSERT(.NOT. ALLOCATED(atom2kind))
207 21 : CPASSERT(.NOT. ALLOCATED(positions))
208 21 : CPASSERT(.NOT. ALLOCATED(xblocks))
209 :
210 21 : natoms = -1
211 21 : nkinds = -1
212 21 : xblocks_read = 0
213 :
214 : CALL open_file(file_name=filename, file_status="OLD", file_form="FORMATTED", &
215 21 : file_action="READ", unit_number=unit_nr)
216 :
217 : ! check if file starts with proper header !TODO: introduce a more unique header
218 21 : READ (unit_nr, fmt=*) label, i1
219 21 : IF (TRIM(label) /= "Version") THEN
220 0 : CPABORT("PAO restart file appears to be corrupted.")
221 : END IF
222 21 : IF (i1 /= file_format_version) CPABORT("Restart PAO file format version is wrong")
223 :
224 : DO WHILE (.TRUE.)
225 377 : READ (unit_nr, fmt=*) label
226 377 : BACKSPACE (unit_nr)
227 :
228 398 : IF (TRIM(label) == "Parametrization") THEN
229 21 : READ (unit_nr, fmt=*) label, str_in
230 21 : param = str_in
231 :
232 356 : ELSE IF (TRIM(label) == "Cell") THEN
233 21 : READ (unit_nr, fmt=*) label, hmat_angstrom
234 21 : ALLOCATE (hmat(3, 3))
235 273 : hmat(:, :) = hmat_angstrom(:, :)/angstrom
236 :
237 335 : ELSE IF (TRIM(label) == "Nkinds") THEN
238 21 : READ (unit_nr, fmt=*) label, nkinds
239 87 : ALLOCATE (kinds(nkinds))
240 :
241 314 : ELSE IF (TRIM(label) == "Kind") THEN
242 24 : READ (unit_nr, fmt=*) label, ikind, str_in, i1
243 24 : CPASSERT(ALLOCATED(kinds))
244 24 : kinds(ikind)%name = str_in
245 24 : kinds(ikind)%z = i1
246 :
247 290 : ELSE IF (TRIM(label) == "PrimBasis") THEN
248 24 : READ (unit_nr, fmt=*) label, ikind, i1, str_in
249 24 : CPASSERT(ALLOCATED(kinds))
250 24 : kinds(ikind)%prim_basis_size = i1
251 24 : kinds(ikind)%prim_basis_name = str_in
252 :
253 266 : ELSE IF (TRIM(label) == "PaoBasis") THEN
254 24 : READ (unit_nr, fmt=*) label, ikind, i1
255 24 : CPASSERT(ALLOCATED(kinds))
256 24 : kinds(ikind)%pao_basis_size = i1
257 :
258 242 : ELSE IF (TRIM(label) == "NPaoPotentials") THEN
259 24 : READ (unit_nr, fmt=*) label, ikind, i1
260 24 : CPASSERT(ALLOCATED(kinds))
261 88 : ALLOCATE (kinds(ikind)%pao_potentials(i1))
262 :
263 218 : ELSE IF (TRIM(label) == "PaoPotential") THEN
264 20 : READ (unit_nr, fmt=*) label, ikind, ipot, i1, i2, r1, r2
265 20 : CPASSERT(ALLOCATED(kinds(ikind)%pao_potentials))
266 20 : kinds(ikind)%pao_potentials(ipot)%maxl = i1
267 20 : kinds(ikind)%pao_potentials(ipot)%max_projector = i2
268 20 : kinds(ikind)%pao_potentials(ipot)%beta = r1
269 20 : kinds(ikind)%pao_potentials(ipot)%weight = r2
270 :
271 198 : ELSE IF (TRIM(label) == "NParams") THEN
272 24 : READ (unit_nr, fmt=*) label, ikind, i1
273 24 : CPASSERT(ALLOCATED(kinds))
274 24 : kinds(ikind)%nparams = i1
275 :
276 174 : ELSE IF (TRIM(label) == "Natoms") THEN
277 21 : READ (unit_nr, fmt=*) label, natoms
278 192 : ALLOCATE (positions(natoms, 3), atom2kind(natoms), xblocks(natoms))
279 66 : positions = 0.0_dp; atom2kind = -1
280 55 : IF (PRESENT(ml_range)) ml_range = [1, natoms]
281 :
282 153 : ELSE IF (TRIM(label) == "MLRange") THEN
283 : ! Natoms entry has to come first
284 0 : CPASSERT(natoms > 0)
285 : ! range of atoms whose xblocks are used for machine learning
286 0 : READ (unit_nr, fmt=*) label, i1, i2
287 0 : IF (PRESENT(ml_range)) ml_range = [i1, i2]
288 :
289 153 : ELSE IF (TRIM(label) == "Atom") THEN
290 45 : READ (unit_nr, fmt=*) label, iatom, str_in, pos_in
291 45 : CPASSERT(ALLOCATED(kinds))
292 51 : DO ikind = 1, nkinds
293 51 : IF (TRIM(kinds(ikind)%name) == TRIM(str_in)) EXIT
294 : END DO
295 45 : CPASSERT(ALLOCATED(atom2kind) .AND. ALLOCATED(positions))
296 45 : atom2kind(iatom) = ikind
297 180 : positions(iatom, :) = pos_in/angstrom
298 :
299 108 : ELSE IF (TRIM(label) == "Xblock") THEN
300 45 : READ (unit_nr, fmt=*) label, iatom
301 45 : CPASSERT(ALLOCATED(kinds) .AND. ALLOCATED(atom2kind))
302 45 : ikind = atom2kind(iatom)
303 45 : nparams = kinds(ikind)%nparams
304 45 : CPASSERT(nparams >= 0)
305 135 : ALLOCATE (xblocks(iatom)%p(nparams, 1))
306 45 : BACKSPACE (unit_nr)
307 499 : READ (unit_nr, fmt=*) label, iatom, xblocks(iatom)%p
308 45 : xblocks_read = xblocks_read + 1
309 45 : CPASSERT(iatom == xblocks_read) ! ensure blocks are read in order
310 :
311 63 : ELSE IF (TRIM(label) == "THE_END") THEN
312 : EXIT
313 : ELSE
314 : !CPWARN("Skipping restart header with label: "//TRIM(label))
315 42 : READ (unit_nr, fmt=*) label ! just read again and ignore
316 : END IF
317 : END DO
318 21 : CALL close_file(unit_number=unit_nr)
319 :
320 21 : CPASSERT(xblocks_read == natoms) ! ensure we read all blocks
321 :
322 21 : END SUBROUTINE pao_read_raw
323 :
324 : ! **************************************************************************************************
325 : !> \brief Ensure that the kind read from the restart is equal to the kind curretly in use.
326 : !> \param pao ...
327 : !> \param qs_env ...
328 : !> \param ikind ...
329 : !> \param pao_kind ...
330 : ! **************************************************************************************************
331 96 : SUBROUTINE pao_kinds_ensure_equal(pao, qs_env, ikind, pao_kind)
332 : TYPE(pao_env_type), POINTER :: pao
333 : TYPE(qs_environment_type), POINTER :: qs_env
334 : INTEGER, INTENT(IN) :: ikind
335 : TYPE(pao_iokind_type), INTENT(IN) :: pao_kind
336 :
337 : CHARACTER(LEN=default_string_length) :: name
338 : INTEGER :: ipot, nparams, pao_basis_size, z
339 24 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
340 : TYPE(gto_basis_set_type), POINTER :: basis_set
341 24 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
342 24 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
343 :
344 : CALL get_qs_env(qs_env, &
345 : atomic_kind_set=atomic_kind_set, &
346 24 : qs_kind_set=qs_kind_set)
347 :
348 24 : IF (ikind > SIZE(atomic_kind_set) .OR. ikind > SIZE(qs_kind_set)) THEN
349 0 : CPABORT("Some kinds are missing.")
350 : END IF
351 :
352 24 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z, name=name)
353 : CALL get_qs_kind(qs_kind_set(ikind), &
354 : basis_set=basis_set, &
355 : pao_basis_size=pao_basis_size, &
356 24 : pao_potentials=pao_potentials)
357 24 : CALL pao_param_count(pao, qs_env, ikind=ikind, nparams=nparams)
358 :
359 24 : IF (pao_kind%nparams /= nparams) THEN
360 0 : CPABORT("Number of parameters do not match")
361 : END IF
362 24 : IF (TRIM(pao_kind%name) /= TRIM(name)) THEN
363 0 : CPABORT("Kind names do not match")
364 : END IF
365 24 : IF (pao_kind%z /= z) THEN
366 0 : CPABORT("Atomic numbers do not match")
367 : END IF
368 24 : IF (TRIM(pao_kind%prim_basis_name) /= TRIM(basis_set%name)) THEN
369 0 : CPABORT("Primary Basis-set name does not match")
370 : END IF
371 24 : IF (pao_kind%prim_basis_size /= basis_set%nsgf) THEN
372 0 : CPABORT("Primary Basis-set size does not match")
373 : END IF
374 24 : IF (pao_kind%pao_basis_size /= pao_basis_size) THEN
375 0 : CPABORT("PAO basis size does not match")
376 : END IF
377 24 : IF (SIZE(pao_kind%pao_potentials) /= SIZE(pao_potentials)) THEN
378 0 : CPABORT("Number of PAO_POTENTIALS does not match")
379 : END IF
380 :
381 44 : DO ipot = 1, SIZE(pao_potentials)
382 20 : IF (pao_kind%pao_potentials(ipot)%maxl /= pao_potentials(ipot)%maxl) THEN
383 0 : CPABORT("PAO_POT_MAXL does not match")
384 : END IF
385 20 : IF (pao_kind%pao_potentials(ipot)%max_projector /= pao_potentials(ipot)%max_projector) THEN
386 0 : CPABORT("PAO_POT_MAX_PROJECTOR does not match")
387 : END IF
388 20 : IF (pao_kind%pao_potentials(ipot)%beta /= pao_potentials(ipot)%beta) THEN
389 0 : CPWARN("PAO_POT_BETA does not match")
390 : END IF
391 44 : IF (pao_kind%pao_potentials(ipot)%weight /= pao_potentials(ipot)%weight) THEN
392 0 : CPWARN("PAO_POT_WEIGHT does not match")
393 : END IF
394 : END DO
395 :
396 24 : END SUBROUTINE pao_kinds_ensure_equal
397 :
398 : ! **************************************************************************************************
399 : !> \brief Writes restart file
400 : !> \param pao ...
401 : !> \param qs_env ...
402 : !> \param energy ...
403 : ! **************************************************************************************************
404 254 : SUBROUTINE pao_write_restart(pao, qs_env, energy)
405 : TYPE(pao_env_type), POINTER :: pao
406 : TYPE(qs_environment_type), POINTER :: qs_env
407 : REAL(dp) :: energy
408 :
409 : CHARACTER(len=*), PARAMETER :: printkey_section = 'DFT%LS_SCF%PAO%PRINT%RESTART', &
410 : routineN = 'pao_write_restart'
411 :
412 : INTEGER :: handle, unit_max, unit_nr
413 : TYPE(cp_logger_type), POINTER :: logger
414 : TYPE(mp_para_env_type), POINTER :: para_env
415 : TYPE(section_vals_type), POINTER :: input
416 :
417 254 : CALL timeset(routineN, handle)
418 254 : logger => cp_get_default_logger()
419 :
420 254 : CALL get_qs_env(qs_env, input=input, para_env=para_env)
421 :
422 : ! open file
423 : unit_nr = cp_print_key_unit_nr(logger, &
424 : input, &
425 : printkey_section, &
426 : extension=".pao", &
427 : file_action="WRITE", &
428 : file_position="REWIND", &
429 : file_status="UNKNOWN", &
430 254 : do_backup=.TRUE.)
431 :
432 : ! although just rank-0 writes the trajectory it requires collective MPI calls
433 254 : unit_max = unit_nr
434 254 : CALL para_env%max(unit_max)
435 254 : IF (unit_max > 0) THEN
436 104 : IF (pao%iw > 0) WRITE (pao%iw, '(A,A)') " PAO| Writing restart file."
437 104 : IF (unit_nr > 0) THEN
438 52 : CALL write_restart_header(pao, qs_env, energy, unit_nr)
439 : END IF
440 :
441 104 : CALL pao_write_diagonal_blocks(para_env, pao%matrix_X, "Xblock", unit_nr)
442 :
443 : END IF
444 :
445 : ! close file
446 254 : IF (unit_nr > 0) WRITE (unit_nr, '(A)') "THE_END"
447 254 : CALL cp_print_key_finished_output(unit_nr, logger, input, printkey_section)
448 :
449 254 : CALL timestop(handle)
450 254 : END SUBROUTINE pao_write_restart
451 :
452 : ! **************************************************************************************************
453 : !> \brief Write the digonal blocks of given DBCSR matrix into the provided unit_nr
454 : !> \param para_env ...
455 : !> \param matrix ...
456 : !> \param label ...
457 : !> \param unit_nr ...
458 : ! **************************************************************************************************
459 104 : SUBROUTINE pao_write_diagonal_blocks(para_env, matrix, label, unit_nr)
460 : TYPE(mp_para_env_type), POINTER :: para_env
461 : TYPE(dbcsr_type) :: matrix
462 : CHARACTER(LEN=*), INTENT(IN) :: label
463 : INTEGER, INTENT(IN) :: unit_nr
464 :
465 : INTEGER :: iatom, natoms
466 104 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
467 : LOGICAL :: found
468 104 : REAL(dp), DIMENSION(:, :), POINTER :: local_block, mpi_buffer
469 :
470 : !TODO: this is a serial algorithm
471 104 : CALL dbcsr_get_info(matrix, row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
472 104 : CPASSERT(SIZE(row_blk_sizes) == SIZE(col_blk_sizes))
473 104 : natoms = SIZE(row_blk_sizes)
474 :
475 352 : DO iatom = 1, natoms
476 984 : ALLOCATE (mpi_buffer(row_blk_sizes(iatom), col_blk_sizes(iatom)))
477 248 : NULLIFY (local_block)
478 248 : CALL dbcsr_get_block_p(matrix=matrix, row=iatom, col=iatom, block=local_block, found=found)
479 248 : IF (ASSOCIATED(local_block)) THEN
480 372 : IF (SIZE(local_block) > 0) THEN
481 : ! catch corner-case
482 4204 : mpi_buffer(:, :) = local_block(:, :)
483 : END IF
484 : ELSE
485 2110 : mpi_buffer(:, :) = 0.0_dp
486 : END IF
487 :
488 8192 : CALL para_env%sum(mpi_buffer)
489 248 : IF (unit_nr > 0) THEN
490 124 : WRITE (unit_nr, fmt="(A,1X,I10,1X)", advance='no') label, iatom
491 2110 : WRITE (unit_nr, *) mpi_buffer
492 : END IF
493 600 : DEALLOCATE (mpi_buffer)
494 : END DO
495 :
496 : ! flush
497 104 : IF (unit_nr > 0) FLUSH (unit_nr)
498 :
499 104 : END SUBROUTINE pao_write_diagonal_blocks
500 :
501 : ! **************************************************************************************************
502 : !> \brief Writes header of restart file
503 : !> \param pao ...
504 : !> \param qs_env ...
505 : !> \param energy ...
506 : !> \param unit_nr ...
507 : ! **************************************************************************************************
508 52 : SUBROUTINE write_restart_header(pao, qs_env, energy, unit_nr)
509 : TYPE(pao_env_type), POINTER :: pao
510 : TYPE(qs_environment_type), POINTER :: qs_env
511 : REAL(dp) :: energy
512 : INTEGER, INTENT(IN) :: unit_nr
513 :
514 : CHARACTER(LEN=default_string_length) :: kindname
515 : INTEGER :: iatom, ikind, ipot, nparams, &
516 : pao_basis_size, z
517 52 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
518 : TYPE(cell_type), POINTER :: cell
519 : TYPE(gto_basis_set_type), POINTER :: basis_set
520 52 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
521 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
522 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
523 :
524 : CALL get_qs_env(qs_env, &
525 : cell=cell, &
526 : particle_set=particle_set, &
527 : atomic_kind_set=atomic_kind_set, &
528 52 : qs_kind_set=qs_kind_set)
529 :
530 52 : WRITE (unit_nr, "(A,5X,I0)") "Version", file_format_version
531 52 : WRITE (unit_nr, "(A,5X,F20.10)") "Energy", energy
532 52 : WRITE (unit_nr, "(A,5X,I0)") "Step", pao%istep
533 52 : WRITE (unit_nr, "(A,5X,A)") "Parametrization", id2str(pao%parameterization)
534 :
535 : ! write kinds
536 52 : WRITE (unit_nr, "(A,5X,I0)") "Nkinds", SIZE(atomic_kind_set)
537 124 : DO ikind = 1, SIZE(atomic_kind_set)
538 72 : CALL get_atomic_kind(atomic_kind_set(ikind), name=kindname, z=z)
539 : CALL get_qs_kind(qs_kind_set(ikind), &
540 : pao_basis_size=pao_basis_size, &
541 : pao_potentials=pao_potentials, &
542 72 : basis_set=basis_set)
543 72 : CALL pao_param_count(pao, qs_env, ikind, nparams)
544 72 : WRITE (unit_nr, "(A,5X,I10,1X,A,1X,I3)") "Kind", ikind, TRIM(kindname), z
545 72 : WRITE (unit_nr, "(A,5X,I10,1X,I3)") "NParams", ikind, nparams
546 72 : WRITE (unit_nr, "(A,5X,I10,1X,I10,1X,A)") "PrimBasis", ikind, basis_set%nsgf, TRIM(basis_set%name)
547 72 : WRITE (unit_nr, "(A,5X,I10,1X,I3)") "PaoBasis", ikind, pao_basis_size
548 72 : WRITE (unit_nr, "(A,5X,I10,1X,I3)") "NPaoPotentials", ikind, SIZE(pao_potentials)
549 245 : DO ipot = 1, SIZE(pao_potentials)
550 49 : WRITE (unit_nr, "(A,5X,I10,1X,I3)", advance='no') "PaoPotential", ikind, ipot
551 49 : WRITE (unit_nr, "(1X,I3)", advance='no') pao_potentials(ipot)%maxl
552 49 : WRITE (unit_nr, "(1X,I3)", advance='no') pao_potentials(ipot)%max_projector
553 49 : WRITE (unit_nr, "(1X,F20.16)", advance='no') pao_potentials(ipot)%beta
554 121 : WRITE (unit_nr, "(1X,F20.16)") pao_potentials(ipot)%weight
555 : END DO
556 : END DO
557 :
558 : ! write cell
559 52 : WRITE (unit_nr, fmt="(A,5X)", advance='no') "Cell"
560 676 : WRITE (unit_nr, *) cell%hmat*angstrom
561 :
562 : ! write atoms
563 52 : WRITE (unit_nr, "(A,5X,I0)") "Natoms", SIZE(particle_set)
564 176 : DO iatom = 1, SIZE(particle_set)
565 124 : kindname = particle_set(iatom)%atomic_kind%name
566 124 : WRITE (unit_nr, fmt="(A,5X,I10,5X,A,1X)", advance='no') "Atom ", iatom, TRIM(kindname)
567 548 : WRITE (unit_nr, *) particle_set(iatom)%r*angstrom
568 : END DO
569 :
570 52 : END SUBROUTINE write_restart_header
571 :
572 : !**************************************************************************************************
573 : !> \brief writing the KS matrix (in terms of the PAO basis) in csr format into a file
574 : !> \param qs_env qs environment
575 : !> \param ls_scf_env ls environment
576 : !> \author Mohammad Hossein Bani-Hashemian
577 : ! **************************************************************************************************
578 294 : SUBROUTINE pao_write_ks_matrix_csr(qs_env, ls_scf_env)
579 : TYPE(qs_environment_type), POINTER :: qs_env
580 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
581 :
582 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_write_ks_matrix_csr'
583 :
584 : CHARACTER(LEN=default_path_length) :: file_name, fileformat
585 : INTEGER :: handle, ispin, output_unit, unit_nr
586 : LOGICAL :: bin, do_kpoints, do_ks_csr_write, uptr
587 : REAL(KIND=dp) :: thld
588 : TYPE(cp_logger_type), POINTER :: logger
589 : TYPE(dbcsr_csr_type) :: ks_mat_csr
590 : TYPE(dbcsr_type) :: matrix_ks_nosym
591 : TYPE(section_vals_type), POINTER :: dft_section, input
592 :
593 294 : CALL timeset(routineN, handle)
594 :
595 294 : NULLIFY (dft_section)
596 :
597 294 : logger => cp_get_default_logger()
598 294 : output_unit = cp_logger_get_default_io_unit(logger)
599 :
600 294 : CALL get_qs_env(qs_env, input=input)
601 294 : dft_section => section_vals_get_subs_vals(input, "DFT")
602 : do_ks_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
603 294 : "PRINT%KS_CSR_WRITE"), cp_p_file)
604 :
605 : ! NOTE: k-points has to be treated differently later. k-points has KS matrix as double pointer.
606 294 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
607 :
608 294 : IF (do_ks_csr_write .AND. (.NOT. do_kpoints)) THEN
609 0 : CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%THRESHOLD", r_val=thld)
610 0 : CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr)
611 0 : CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%BINARY", l_val=bin)
612 :
613 0 : IF (bin) THEN
614 0 : fileformat = "UNFORMATTED"
615 : ELSE
616 0 : fileformat = "FORMATTED"
617 : END IF
618 :
619 0 : DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
620 :
621 0 : IF (dbcsr_has_symmetry(ls_scf_env%matrix_ks(ispin))) THEN
622 0 : CALL dbcsr_desymmetrize(ls_scf_env%matrix_ks(ispin), matrix_ks_nosym)
623 : ELSE
624 0 : CALL dbcsr_copy(matrix_ks_nosym, ls_scf_env%matrix_ks(ispin))
625 : END IF
626 :
627 0 : CALL dbcsr_csr_create_from_dbcsr(matrix_ks_nosym, ks_mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
628 0 : CALL dbcsr_convert_dbcsr_to_csr(matrix_ks_nosym, ks_mat_csr)
629 :
630 0 : WRITE (file_name, '(A,I0)') "PAO_KS_SPIN_", ispin
631 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%KS_CSR_WRITE", &
632 : extension=".csr", middle_name=TRIM(file_name), &
633 0 : file_status="REPLACE", file_form=fileformat)
634 0 : CALL dbcsr_csr_write(ks_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
635 :
636 0 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%KS_CSR_WRITE")
637 :
638 0 : CALL dbcsr_csr_destroy(ks_mat_csr)
639 0 : CALL dbcsr_release(matrix_ks_nosym)
640 : END DO
641 : END IF
642 :
643 294 : CALL timestop(handle)
644 :
645 294 : END SUBROUTINE pao_write_ks_matrix_csr
646 :
647 : !**************************************************************************************************
648 : !> \brief writing the overlap matrix (in terms of the PAO basis) in csr format into a file
649 : !> \param qs_env qs environment
650 : !> \param ls_scf_env ls environment
651 : !> \author Mohammad Hossein Bani-Hashemian
652 : ! **************************************************************************************************
653 294 : SUBROUTINE pao_write_s_matrix_csr(qs_env, ls_scf_env)
654 : TYPE(qs_environment_type), POINTER :: qs_env
655 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
656 :
657 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_write_s_matrix_csr'
658 :
659 : CHARACTER(LEN=default_path_length) :: file_name, fileformat
660 : INTEGER :: handle, output_unit, unit_nr
661 : LOGICAL :: bin, do_kpoints, do_s_csr_write, uptr
662 : REAL(KIND=dp) :: thld
663 : TYPE(cp_logger_type), POINTER :: logger
664 : TYPE(dbcsr_csr_type) :: s_mat_csr
665 : TYPE(dbcsr_type) :: matrix_s_nosym
666 : TYPE(section_vals_type), POINTER :: dft_section, input
667 :
668 294 : CALL timeset(routineN, handle)
669 :
670 294 : NULLIFY (dft_section)
671 :
672 294 : logger => cp_get_default_logger()
673 294 : output_unit = cp_logger_get_default_io_unit(logger)
674 :
675 294 : CALL get_qs_env(qs_env, input=input)
676 294 : dft_section => section_vals_get_subs_vals(input, "DFT")
677 : do_s_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
678 294 : "PRINT%S_CSR_WRITE"), cp_p_file)
679 :
680 : ! NOTE: k-points has to be treated differently later. k-points has overlap matrix as double pointer.
681 294 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
682 :
683 294 : IF (do_s_csr_write .AND. (.NOT. do_kpoints)) THEN
684 0 : CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%THRESHOLD", r_val=thld)
685 0 : CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr)
686 0 : CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%BINARY", l_val=bin)
687 :
688 0 : IF (bin) THEN
689 0 : fileformat = "UNFORMATTED"
690 : ELSE
691 0 : fileformat = "FORMATTED"
692 : END IF
693 :
694 0 : IF (dbcsr_has_symmetry(ls_scf_env%matrix_s)) THEN
695 0 : CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_s_nosym)
696 : ELSE
697 0 : CALL dbcsr_copy(matrix_s_nosym, ls_scf_env%matrix_s)
698 : END IF
699 :
700 0 : CALL dbcsr_csr_create_from_dbcsr(matrix_s_nosym, s_mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
701 0 : CALL dbcsr_convert_dbcsr_to_csr(matrix_s_nosym, s_mat_csr)
702 :
703 0 : WRITE (file_name, '(A,I0)') "PAO_S"
704 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%S_CSR_WRITE", &
705 : extension=".csr", middle_name=TRIM(file_name), &
706 0 : file_status="REPLACE", file_form=fileformat)
707 0 : CALL dbcsr_csr_write(s_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
708 :
709 0 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%S_CSR_WRITE")
710 :
711 0 : CALL dbcsr_csr_destroy(s_mat_csr)
712 0 : CALL dbcsr_release(matrix_s_nosym)
713 : END IF
714 :
715 294 : CALL timestop(handle)
716 :
717 294 : END SUBROUTINE pao_write_s_matrix_csr
718 :
719 : !**************************************************************************************************
720 : !> \brief writing the core Hamiltonian matrix (NYA)
721 : !> \param qs_env qs environment
722 : !> \param ls_scf_env ls environment
723 : !> \author Mohammad Hossein Bani-Hashemian
724 : ! **************************************************************************************************
725 294 : SUBROUTINE pao_write_hcore_matrix_csr(qs_env, ls_scf_env)
726 : TYPE(qs_environment_type), POINTER :: qs_env
727 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
728 :
729 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_write_hcore_matrix_csr'
730 :
731 : INTEGER :: handle, output_unit
732 : LOGICAL :: do_h_csr_write, do_kpoints
733 : TYPE(cp_logger_type), POINTER :: logger
734 : TYPE(section_vals_type), POINTER :: dft_section, input
735 :
736 : MARK_USED(ls_scf_env)
737 :
738 294 : CALL timeset(routineN, handle)
739 :
740 294 : NULLIFY (dft_section)
741 :
742 294 : logger => cp_get_default_logger()
743 294 : output_unit = cp_logger_get_default_io_unit(logger)
744 :
745 294 : CALL get_qs_env(qs_env, input=input)
746 294 : dft_section => section_vals_get_subs_vals(input, "DFT")
747 : do_h_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
748 294 : "PRINT%HCORE_CSR_WRITE"), cp_p_file)
749 :
750 : ! NOTE: k-points has to be treated differently later. k-points has KS matrix as double pointer.
751 294 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
752 :
753 294 : IF (do_h_csr_write .AND. (.NOT. do_kpoints)) THEN
754 0 : CALL cp_warn(__LOCATION__, "Writing the PAO Core Hamiltonian matrix in CSR format NYA")
755 : END IF
756 :
757 294 : CALL timestop(handle)
758 :
759 294 : END SUBROUTINE pao_write_hcore_matrix_csr
760 :
761 : !**************************************************************************************************
762 : !> \brief writing the density matrix (NYA)
763 : !> \param qs_env qs environment
764 : !> \param ls_scf_env ls environment
765 : !> \author Mohammad Hossein Bani-Hashemian
766 : ! **************************************************************************************************
767 294 : SUBROUTINE pao_write_p_matrix_csr(qs_env, ls_scf_env)
768 : TYPE(qs_environment_type), POINTER :: qs_env
769 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
770 :
771 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_write_p_matrix_csr'
772 :
773 : INTEGER :: handle, output_unit
774 : LOGICAL :: do_kpoints, do_p_csr_write
775 : TYPE(cp_logger_type), POINTER :: logger
776 : TYPE(section_vals_type), POINTER :: dft_section, input
777 :
778 : MARK_USED(ls_scf_env)
779 :
780 294 : CALL timeset(routineN, handle)
781 :
782 294 : NULLIFY (dft_section)
783 :
784 294 : logger => cp_get_default_logger()
785 294 : output_unit = cp_logger_get_default_io_unit(logger)
786 :
787 294 : CALL get_qs_env(qs_env, input=input)
788 294 : dft_section => section_vals_get_subs_vals(input, "DFT")
789 : do_p_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
790 294 : "PRINT%P_CSR_WRITE"), cp_p_file)
791 :
792 : ! NOTE: k-points has to be treated differently later. k-points has KS matrix as double pointer.
793 294 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
794 :
795 294 : IF (do_p_csr_write .AND. (.NOT. do_kpoints)) THEN
796 0 : CALL cp_warn(__LOCATION__, "Writing the PAO density matrix in CSR format NYA")
797 : END IF
798 :
799 294 : CALL timestop(handle)
800 :
801 294 : END SUBROUTINE pao_write_p_matrix_csr
802 :
803 0 : END MODULE pao_io
|
