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