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