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