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_files, ONLY: close_file,&
18 : open_file
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_get_default_io_unit,&
21 : cp_logger_type
22 : USE cp_output_handling, ONLY: cp_p_file,&
23 : cp_print_key_finished_output,&
24 : cp_print_key_should_output,&
25 : cp_print_key_unit_nr
26 : USE dbcsr_api, ONLY: &
27 : dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_csr_create_from_dbcsr, &
28 : dbcsr_csr_dbcsr_blkrow_dist, dbcsr_csr_destroy, dbcsr_csr_type, dbcsr_csr_write, &
29 : dbcsr_desymmetrize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_has_symmetry, dbcsr_release, &
30 : dbcsr_type
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) &
119 0 : CPWARN("Restarting from different cell")
120 :
121 : ! check parametrization
122 4 : IF (TRIM(param) .NE. TRIM(ADJUSTL(id2str(pao%parameterization)))) &
123 4 : CPABORT("Restart PAO parametrization does not match")
124 :
125 : ! check kinds
126 11 : DO ikind = 1, SIZE(kinds)
127 11 : CALL pao_kinds_ensure_equal(pao, qs_env, ikind, kinds(ikind))
128 : END DO
129 :
130 : ! check number of atoms
131 4 : IF (SIZE(positions, 1) /= natoms) &
132 0 : CPABORT("Number of atoms do not match")
133 :
134 : ! check atom2kind
135 15 : DO iatom = 1, natoms
136 11 : IF (atom2kind(iatom) /= particle_set(iatom)%atomic_kind%kind_number) &
137 4 : CPABORT("Restart atomic kinds do not match.")
138 : END DO
139 :
140 : ! check positions, warning only
141 4 : diff = 0.0_dp
142 15 : DO iatom = 1, natoms
143 59 : diff = MAX(diff, MAXVAL(ABS(positions(iatom, :) - particle_set(iatom)%r)))
144 : END DO
145 4 : IF (diff > 1e-10) &
146 1 : CPWARN("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) &
369 0 : CPABORT("PAO_POT_MAXL does not match")
370 20 : IF (pao_kind%pao_potentials(ipot)%max_projector /= pao_potentials(ipot)%max_projector) &
371 0 : CPABORT("PAO_POT_MAX_PROJECTOR does not match")
372 20 : IF (pao_kind%pao_potentials(ipot)%beta /= pao_potentials(ipot)%beta) &
373 0 : CPWARN("PAO_POT_BETA does not match")
374 20 : IF (pao_kind%pao_potentials(ipot)%weight /= pao_potentials(ipot)%weight) &
375 24 : CPWARN("PAO_POT_WEIGHT does not match")
376 : END DO
377 :
378 24 : END SUBROUTINE pao_kinds_ensure_equal
379 :
380 : ! **************************************************************************************************
381 : !> \brief Writes restart file
382 : !> \param pao ...
383 : !> \param qs_env ...
384 : !> \param energy ...
385 : ! **************************************************************************************************
386 254 : SUBROUTINE pao_write_restart(pao, qs_env, energy)
387 : TYPE(pao_env_type), POINTER :: pao
388 : TYPE(qs_environment_type), POINTER :: qs_env
389 : REAL(dp) :: energy
390 :
391 : CHARACTER(len=*), PARAMETER :: printkey_section = 'DFT%LS_SCF%PAO%PRINT%RESTART', &
392 : routineN = 'pao_write_restart'
393 :
394 : INTEGER :: handle, unit_max, unit_nr
395 : TYPE(cp_logger_type), POINTER :: logger
396 : TYPE(mp_para_env_type), POINTER :: para_env
397 : TYPE(section_vals_type), POINTER :: input
398 :
399 254 : CALL timeset(routineN, handle)
400 254 : logger => cp_get_default_logger()
401 :
402 254 : CALL get_qs_env(qs_env, input=input, para_env=para_env)
403 :
404 : ! open file
405 : unit_nr = cp_print_key_unit_nr(logger, &
406 : input, &
407 : printkey_section, &
408 : extension=".pao", &
409 : file_action="WRITE", &
410 : file_position="REWIND", &
411 : file_status="UNKNOWN", &
412 254 : do_backup=.TRUE.)
413 :
414 : ! although just rank-0 writes the trajectory it requires collective MPI calls
415 254 : unit_max = unit_nr
416 254 : CALL para_env%max(unit_max)
417 254 : IF (unit_max > 0) THEN
418 104 : IF (pao%iw > 0) WRITE (pao%iw, '(A,A)') " PAO| Writing restart file."
419 104 : IF (unit_nr > 0) &
420 52 : CALL write_restart_header(pao, qs_env, energy, unit_nr)
421 :
422 104 : CALL pao_write_diagonal_blocks(para_env, pao%matrix_X, "Xblock", unit_nr)
423 :
424 : END IF
425 :
426 : ! close file
427 254 : IF (unit_nr > 0) WRITE (unit_nr, '(A)') "THE_END"
428 254 : CALL cp_print_key_finished_output(unit_nr, logger, input, printkey_section)
429 :
430 254 : CALL timestop(handle)
431 254 : END SUBROUTINE pao_write_restart
432 :
433 : ! **************************************************************************************************
434 : !> \brief Write the digonal blocks of given DBCSR matrix into the provided unit_nr
435 : !> \param para_env ...
436 : !> \param matrix ...
437 : !> \param label ...
438 : !> \param unit_nr ...
439 : ! **************************************************************************************************
440 104 : SUBROUTINE pao_write_diagonal_blocks(para_env, matrix, label, unit_nr)
441 : TYPE(mp_para_env_type), POINTER :: para_env
442 : TYPE(dbcsr_type) :: matrix
443 : CHARACTER(LEN=*), INTENT(IN) :: label
444 : INTEGER, INTENT(IN) :: unit_nr
445 :
446 : INTEGER :: iatom, natoms
447 104 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
448 : LOGICAL :: found
449 104 : REAL(dp), DIMENSION(:, :), POINTER :: local_block, mpi_buffer
450 :
451 : !TODO: this is a serial algorithm
452 104 : CALL dbcsr_get_info(matrix, row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
453 104 : CPASSERT(SIZE(row_blk_sizes) == SIZE(col_blk_sizes))
454 104 : natoms = SIZE(row_blk_sizes)
455 :
456 352 : DO iatom = 1, natoms
457 984 : ALLOCATE (mpi_buffer(row_blk_sizes(iatom), col_blk_sizes(iatom)))
458 248 : NULLIFY (local_block)
459 248 : CALL dbcsr_get_block_p(matrix=matrix, row=iatom, col=iatom, block=local_block, found=found)
460 248 : IF (ASSOCIATED(local_block)) THEN
461 124 : IF (SIZE(local_block) > 0) & ! catch corner-case
462 4084 : mpi_buffer(:, :) = local_block(:, :)
463 : ELSE
464 2110 : mpi_buffer(:, :) = 0.0_dp
465 : END IF
466 :
467 8192 : CALL para_env%sum(mpi_buffer)
468 248 : IF (unit_nr > 0) THEN
469 124 : WRITE (unit_nr, fmt="(A,1X,I10,1X)", advance='no') label, iatom
470 2110 : WRITE (unit_nr, *) mpi_buffer
471 : END IF
472 600 : DEALLOCATE (mpi_buffer)
473 : END DO
474 :
475 : ! flush
476 104 : IF (unit_nr > 0) FLUSH (unit_nr)
477 :
478 104 : END SUBROUTINE pao_write_diagonal_blocks
479 :
480 : ! **************************************************************************************************
481 : !> \brief Writes header of restart file
482 : !> \param pao ...
483 : !> \param qs_env ...
484 : !> \param energy ...
485 : !> \param unit_nr ...
486 : ! **************************************************************************************************
487 52 : SUBROUTINE write_restart_header(pao, qs_env, energy, unit_nr)
488 : TYPE(pao_env_type), POINTER :: pao
489 : TYPE(qs_environment_type), POINTER :: qs_env
490 : REAL(dp) :: energy
491 : INTEGER, INTENT(IN) :: unit_nr
492 :
493 : CHARACTER(LEN=default_string_length) :: kindname
494 : INTEGER :: iatom, ikind, ipot, nparams, &
495 : pao_basis_size, z
496 52 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
497 : TYPE(cell_type), POINTER :: cell
498 : TYPE(gto_basis_set_type), POINTER :: basis_set
499 52 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
500 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
501 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
502 :
503 : CALL get_qs_env(qs_env, &
504 : cell=cell, &
505 : particle_set=particle_set, &
506 : atomic_kind_set=atomic_kind_set, &
507 52 : qs_kind_set=qs_kind_set)
508 :
509 52 : WRITE (unit_nr, "(A,5X,I0)") "Version", file_format_version
510 52 : WRITE (unit_nr, "(A,5X,F20.10)") "Energy", energy
511 52 : WRITE (unit_nr, "(A,5X,I0)") "Step", pao%istep
512 52 : WRITE (unit_nr, "(A,5X,A)") "Parametrization", id2str(pao%parameterization)
513 :
514 : ! write kinds
515 52 : WRITE (unit_nr, "(A,5X,I0)") "Nkinds", SIZE(atomic_kind_set)
516 124 : DO ikind = 1, SIZE(atomic_kind_set)
517 72 : CALL get_atomic_kind(atomic_kind_set(ikind), name=kindname, z=z)
518 : CALL get_qs_kind(qs_kind_set(ikind), &
519 : pao_basis_size=pao_basis_size, &
520 : pao_potentials=pao_potentials, &
521 72 : basis_set=basis_set)
522 72 : CALL pao_param_count(pao, qs_env, ikind, nparams)
523 72 : WRITE (unit_nr, "(A,5X,I10,1X,A,1X,I3)") "Kind", ikind, TRIM(kindname), z
524 72 : WRITE (unit_nr, "(A,5X,I10,1X,I3)") "NParams", ikind, nparams
525 72 : WRITE (unit_nr, "(A,5X,I10,1X,I10,1X,A)") "PrimBasis", ikind, basis_set%nsgf, TRIM(basis_set%name)
526 72 : WRITE (unit_nr, "(A,5X,I10,1X,I3)") "PaoBasis", ikind, pao_basis_size
527 72 : WRITE (unit_nr, "(A,5X,I10,1X,I3)") "NPaoPotentials", ikind, SIZE(pao_potentials)
528 245 : DO ipot = 1, SIZE(pao_potentials)
529 49 : WRITE (unit_nr, "(A,5X,I10,1X,I3)", advance='no') "PaoPotential", ikind, ipot
530 49 : WRITE (unit_nr, "(1X,I3)", advance='no') pao_potentials(ipot)%maxl
531 49 : WRITE (unit_nr, "(1X,I3)", advance='no') pao_potentials(ipot)%max_projector
532 49 : WRITE (unit_nr, "(1X,F20.16)", advance='no') pao_potentials(ipot)%beta
533 121 : WRITE (unit_nr, "(1X,F20.16)") pao_potentials(ipot)%weight
534 : END DO
535 : END DO
536 :
537 : ! write cell
538 52 : WRITE (unit_nr, fmt="(A,5X)", advance='no') "Cell"
539 676 : WRITE (unit_nr, *) cell%hmat*angstrom
540 :
541 : ! write atoms
542 52 : WRITE (unit_nr, "(A,5X,I0)") "Natoms", SIZE(particle_set)
543 176 : DO iatom = 1, SIZE(particle_set)
544 124 : kindname = particle_set(iatom)%atomic_kind%name
545 124 : WRITE (unit_nr, fmt="(A,5X,I10,5X,A,1X)", advance='no') "Atom ", iatom, TRIM(kindname)
546 548 : WRITE (unit_nr, *) particle_set(iatom)%r*angstrom
547 : END DO
548 :
549 52 : END SUBROUTINE write_restart_header
550 :
551 : !**************************************************************************************************
552 : !> \brief writing the KS matrix (in terms of the PAO basis) in csr format into a file
553 : !> \param qs_env qs environment
554 : !> \param ls_scf_env ls environment
555 : !> \author Mohammad Hossein Bani-Hashemian
556 : ! **************************************************************************************************
557 278 : SUBROUTINE pao_write_ks_matrix_csr(qs_env, ls_scf_env)
558 : TYPE(qs_environment_type), POINTER :: qs_env
559 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
560 :
561 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_write_ks_matrix_csr'
562 :
563 : CHARACTER(LEN=default_path_length) :: file_name, fileformat
564 : INTEGER :: handle, ispin, output_unit, unit_nr
565 : LOGICAL :: bin, do_kpoints, do_ks_csr_write, uptr
566 : REAL(KIND=dp) :: thld
567 : TYPE(cp_logger_type), POINTER :: logger
568 : TYPE(dbcsr_csr_type) :: ks_mat_csr
569 : TYPE(dbcsr_type) :: matrix_ks_nosym
570 : TYPE(section_vals_type), POINTER :: dft_section, input
571 :
572 278 : CALL timeset(routineN, handle)
573 :
574 278 : NULLIFY (dft_section)
575 :
576 278 : logger => cp_get_default_logger()
577 278 : output_unit = cp_logger_get_default_io_unit(logger)
578 :
579 278 : CALL get_qs_env(qs_env, input=input)
580 278 : dft_section => section_vals_get_subs_vals(input, "DFT")
581 : do_ks_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
582 278 : "PRINT%KS_CSR_WRITE"), cp_p_file)
583 :
584 : ! NOTE: k-points has to be treated differently later. k-points has KS matrix as double pointer.
585 278 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
586 :
587 278 : IF (do_ks_csr_write .AND. (.NOT. do_kpoints)) THEN
588 0 : CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%THRESHOLD", r_val=thld)
589 0 : CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr)
590 0 : CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%BINARY", l_val=bin)
591 :
592 0 : IF (bin) THEN
593 0 : fileformat = "UNFORMATTED"
594 : ELSE
595 0 : fileformat = "FORMATTED"
596 : END IF
597 :
598 0 : DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
599 :
600 0 : IF (dbcsr_has_symmetry(ls_scf_env%matrix_ks(ispin))) THEN
601 0 : CALL dbcsr_desymmetrize(ls_scf_env%matrix_ks(ispin), matrix_ks_nosym)
602 : ELSE
603 0 : CALL dbcsr_copy(matrix_ks_nosym, ls_scf_env%matrix_ks(ispin))
604 : END IF
605 :
606 0 : CALL dbcsr_csr_create_from_dbcsr(matrix_ks_nosym, ks_mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
607 0 : CALL dbcsr_convert_dbcsr_to_csr(matrix_ks_nosym, ks_mat_csr)
608 :
609 0 : WRITE (file_name, '(A,I0)') "PAO_KS_SPIN_", ispin
610 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%KS_CSR_WRITE", &
611 : extension=".csr", middle_name=TRIM(file_name), &
612 0 : file_status="REPLACE", file_form=fileformat)
613 0 : CALL dbcsr_csr_write(ks_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
614 :
615 0 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%KS_CSR_WRITE")
616 :
617 0 : CALL dbcsr_csr_destroy(ks_mat_csr)
618 0 : CALL dbcsr_release(matrix_ks_nosym)
619 : END DO
620 : END IF
621 :
622 278 : CALL timestop(handle)
623 :
624 278 : END SUBROUTINE pao_write_ks_matrix_csr
625 :
626 : !**************************************************************************************************
627 : !> \brief writing the overlap matrix (in terms of the PAO basis) in csr format into a file
628 : !> \param qs_env qs environment
629 : !> \param ls_scf_env ls environment
630 : !> \author Mohammad Hossein Bani-Hashemian
631 : ! **************************************************************************************************
632 278 : SUBROUTINE pao_write_s_matrix_csr(qs_env, ls_scf_env)
633 : TYPE(qs_environment_type), POINTER :: qs_env
634 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
635 :
636 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_write_s_matrix_csr'
637 :
638 : CHARACTER(LEN=default_path_length) :: file_name, fileformat
639 : INTEGER :: handle, output_unit, unit_nr
640 : LOGICAL :: bin, do_kpoints, do_s_csr_write, uptr
641 : REAL(KIND=dp) :: thld
642 : TYPE(cp_logger_type), POINTER :: logger
643 : TYPE(dbcsr_csr_type) :: s_mat_csr
644 : TYPE(dbcsr_type) :: matrix_s_nosym
645 : TYPE(section_vals_type), POINTER :: dft_section, input
646 :
647 278 : CALL timeset(routineN, handle)
648 :
649 278 : NULLIFY (dft_section)
650 :
651 278 : logger => cp_get_default_logger()
652 278 : output_unit = cp_logger_get_default_io_unit(logger)
653 :
654 278 : CALL get_qs_env(qs_env, input=input)
655 278 : dft_section => section_vals_get_subs_vals(input, "DFT")
656 : do_s_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
657 278 : "PRINT%S_CSR_WRITE"), cp_p_file)
658 :
659 : ! NOTE: k-points has to be treated differently later. k-points has overlap matrix as double pointer.
660 278 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
661 :
662 278 : IF (do_s_csr_write .AND. (.NOT. do_kpoints)) THEN
663 0 : CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%THRESHOLD", r_val=thld)
664 0 : CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr)
665 0 : CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%BINARY", l_val=bin)
666 :
667 0 : IF (bin) THEN
668 0 : fileformat = "UNFORMATTED"
669 : ELSE
670 0 : fileformat = "FORMATTED"
671 : END IF
672 :
673 0 : IF (dbcsr_has_symmetry(ls_scf_env%matrix_s)) THEN
674 0 : CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_s_nosym)
675 : ELSE
676 0 : CALL dbcsr_copy(matrix_s_nosym, ls_scf_env%matrix_s)
677 : END IF
678 :
679 0 : CALL dbcsr_csr_create_from_dbcsr(matrix_s_nosym, s_mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
680 0 : CALL dbcsr_convert_dbcsr_to_csr(matrix_s_nosym, s_mat_csr)
681 :
682 0 : WRITE (file_name, '(A,I0)') "PAO_S"
683 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%S_CSR_WRITE", &
684 : extension=".csr", middle_name=TRIM(file_name), &
685 0 : file_status="REPLACE", file_form=fileformat)
686 0 : CALL dbcsr_csr_write(s_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
687 :
688 0 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%S_CSR_WRITE")
689 :
690 0 : CALL dbcsr_csr_destroy(s_mat_csr)
691 0 : CALL dbcsr_release(matrix_s_nosym)
692 : END IF
693 :
694 278 : CALL timestop(handle)
695 :
696 278 : END SUBROUTINE pao_write_s_matrix_csr
697 :
698 0 : END MODULE pao_io
|
