(git:374b731)
Loading...
Searching...
No Matches
pao_io.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE pao_io
16 USE cell_types, ONLY: cell_type
17 USE cp_files, ONLY: close_file,&
22 USE cp_output_handling, ONLY: cp_p_file,&
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
35 USE kinds, ONLY: default_path_length,&
37 dp
39 USE pao_input, ONLY: id2str
40 USE pao_param, ONLY: pao_param_count
41 USE pao_types, ONLY: pao_env_type
43 USE physcon, ONLY: angstrom
46 USE qs_kind_types, ONLY: get_qs_kind,&
49#include "./base/base_uses.f90"
50
51 IMPLICIT NONE
52
53 PRIVATE
54
55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_io'
56
61
62 ! data types used by pao_read_raw()
64 REAL(dp), DIMENSION(:, :), ALLOCATABLE :: p
65 END TYPE pao_ioblock_type
66
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
79CONTAINS
80
81! **************************************************************************************************
82!> \brief Reads restart file
83!> \param pao ...
84!> \param qs_env ...
85! **************************************************************************************************
86 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2kind
93 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
94 LOGICAL :: found
95 REAL(dp) :: diff
96 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hmat, positions
97 REAL(dp), DIMENSION(:, :), POINTER :: block_x, buffer
98 TYPE(cell_type), POINTER :: cell
99 TYPE(mp_para_env_type), POINTER :: para_env
100 TYPE(pao_ioblock_type), ALLOCATABLE, DIMENSION(:) :: xblocks
101 TYPE(pao_iokind_type), ALLOCATABLE, DIMENSION(:) :: kinds
102 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
103
104 cpassert(len_trim(pao%restart_file) > 0)
105 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 particle_set=particle_set)
112
113 ! read and check restart file on first rank only
114 IF (para_env%is_source()) THEN
115 CALL pao_read_raw(pao%restart_file, param, hmat, kinds, atom2kind, positions, xblocks)
116
117 ! check cell
118 IF (maxval(abs(hmat - cell%hmat)) > 1e-10) &
119 cpwarn("Restarting from different cell")
120
121 ! check parametrization
122 IF (trim(param) .NE. trim(adjustl(id2str(pao%parameterization)))) &
123 cpabort("Restart PAO parametrization does not match")
124
125 ! check kinds
126 DO ikind = 1, SIZE(kinds)
127 CALL pao_kinds_ensure_equal(pao, qs_env, ikind, kinds(ikind))
128 END DO
129
130 ! check number of atoms
131 IF (SIZE(positions, 1) /= natoms) &
132 cpabort("Number of atoms do not match")
133
134 ! check atom2kind
135 DO iatom = 1, natoms
136 IF (atom2kind(iatom) /= particle_set(iatom)%atomic_kind%kind_number) &
137 cpabort("Restart atomic kinds do not match.")
138 END DO
139
140 ! check positions, warning only
141 diff = 0.0_dp
142 DO iatom = 1, natoms
143 diff = max(diff, maxval(abs(positions(iatom, :) - particle_set(iatom)%r)))
144 END DO
145 IF (diff > 1e-10) &
146 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 CALL dbcsr_get_info(pao%matrix_X, row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
153 DO iatom = 1, natoms
154 ALLOCATE (buffer(row_blk_sizes(iatom), col_blk_sizes(iatom)))
155 IF (para_env%is_source()) THEN
156 cpassert(row_blk_sizes(iatom) == SIZE(xblocks(iatom)%p, 1))
157 cpassert(col_blk_sizes(iatom) == SIZE(xblocks(iatom)%p, 2))
158 buffer = xblocks(iatom)%p
159 END IF
160 CALL para_env%bcast(buffer)
161 CALL dbcsr_get_block_p(matrix=pao%matrix_X, row=iatom, col=iatom, block=block_x, found=found)
162 IF (ASSOCIATED(block_x)) &
163 block_x = buffer
164 DEALLOCATE (buffer)
165 END DO
166
167 ! ALLOCATABLEs deallocate themselves
168
169 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 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 cpassert(.NOT. ALLOCATED(hmat))
200 cpassert(.NOT. ALLOCATED(kinds))
201 cpassert(.NOT. ALLOCATED(atom2kind))
202 cpassert(.NOT. ALLOCATED(positions))
203 cpassert(.NOT. ALLOCATED(xblocks))
204
205 natoms = -1
206 nkinds = -1
207 xblocks_read = 0
208
209 CALL open_file(file_name=filename, file_status="OLD", file_form="FORMATTED", &
210 file_action="READ", unit_number=unit_nr)
211
212 ! check if file starts with proper header !TODO: introduce a more unique header
213 READ (unit_nr, fmt=*) label, i1
214 IF (trim(label) /= "Version") &
215 cpabort("PAO restart file appears to be corrupted.")
216 IF (i1 /= file_format_version) cpabort("Restart PAO file format version is wrong")
217
218 DO WHILE (.true.)
219 READ (unit_nr, fmt=*) label
220 backspace(unit_nr)
221
222 IF (trim(label) == "Parametrization") THEN
223 READ (unit_nr, fmt=*) label, str_in
224 param = str_in
225
226 ELSE IF (trim(label) == "Cell") THEN
227 READ (unit_nr, fmt=*) label, hmat_angstrom
228 ALLOCATE (hmat(3, 3))
229 hmat(:, :) = hmat_angstrom(:, :)/angstrom
230
231 ELSE IF (trim(label) == "Nkinds") THEN
232 READ (unit_nr, fmt=*) label, nkinds
233 ALLOCATE (kinds(nkinds))
234
235 ELSE IF (trim(label) == "Kind") THEN
236 READ (unit_nr, fmt=*) label, ikind, str_in, i1
237 cpassert(ALLOCATED(kinds))
238 kinds(ikind)%name = str_in
239 kinds(ikind)%z = i1
240
241 ELSE IF (trim(label) == "PrimBasis") THEN
242 READ (unit_nr, fmt=*) label, ikind, i1, str_in
243 cpassert(ALLOCATED(kinds))
244 kinds(ikind)%prim_basis_size = i1
245 kinds(ikind)%prim_basis_name = str_in
246
247 ELSE IF (trim(label) == "PaoBasis") THEN
248 READ (unit_nr, fmt=*) label, ikind, i1
249 cpassert(ALLOCATED(kinds))
250 kinds(ikind)%pao_basis_size = i1
251
252 ELSE IF (trim(label) == "NPaoPotentials") THEN
253 READ (unit_nr, fmt=*) label, ikind, i1
254 cpassert(ALLOCATED(kinds))
255 ALLOCATE (kinds(ikind)%pao_potentials(i1))
256
257 ELSE IF (trim(label) == "PaoPotential") THEN
258 READ (unit_nr, fmt=*) label, ikind, ipot, i1, i2, r1, r2
259 cpassert(ALLOCATED(kinds(ikind)%pao_potentials))
260 kinds(ikind)%pao_potentials(ipot)%maxl = i1
261 kinds(ikind)%pao_potentials(ipot)%max_projector = i2
262 kinds(ikind)%pao_potentials(ipot)%beta = r1
263 kinds(ikind)%pao_potentials(ipot)%weight = r2
264
265 ELSE IF (trim(label) == "NParams") THEN
266 READ (unit_nr, fmt=*) label, ikind, i1
267 cpassert(ALLOCATED(kinds))
268 kinds(ikind)%nparams = i1
269
270 ELSE IF (trim(label) == "Natoms") THEN
271 READ (unit_nr, fmt=*) label, natoms
272 ALLOCATE (positions(natoms, 3), atom2kind(natoms), xblocks(natoms))
273 positions = 0.0_dp; atom2kind = -1
274 IF (PRESENT(ml_range)) ml_range = (/1, natoms/)
275
276 ELSE IF (trim(label) == "MLRange") THEN
277 ! Natoms entry has to come first
278 cpassert(natoms > 0)
279 ! range of atoms whose xblocks are used for machine learning
280 READ (unit_nr, fmt=*) label, i1, i2
281 IF (PRESENT(ml_range)) ml_range = (/i1, i2/)
282
283 ELSE IF (trim(label) == "Atom") THEN
284 READ (unit_nr, fmt=*) label, iatom, str_in, pos_in
285 cpassert(ALLOCATED(kinds))
286 DO ikind = 1, nkinds
287 IF (trim(kinds(ikind)%name) .EQ. trim(str_in)) EXIT
288 END DO
289 cpassert(ALLOCATED(atom2kind) .AND. ALLOCATED(positions))
290 atom2kind(iatom) = ikind
291 positions(iatom, :) = pos_in/angstrom
292
293 ELSE IF (trim(label) == "Xblock") THEN
294 READ (unit_nr, fmt=*) label, iatom
295 cpassert(ALLOCATED(kinds) .AND. ALLOCATED(atom2kind))
296 ikind = atom2kind(iatom)
297 nparams = kinds(ikind)%nparams
298 cpassert(nparams >= 0)
299 ALLOCATE (xblocks(iatom)%p(nparams, 1))
300 backspace(unit_nr)
301 READ (unit_nr, fmt=*) label, iatom, xblocks(iatom)%p
302 xblocks_read = xblocks_read + 1
303 cpassert(iatom == xblocks_read) ! ensure blocks are read in order
304
305 ELSE IF (trim(label) == "THE_END") THEN
306 EXIT
307 ELSE
308 !CPWARN("Skipping restart header with label: "//TRIM(label))
309 READ (unit_nr, fmt=*) label ! just read again and ignore
310 END IF
311 END DO
312 CALL close_file(unit_number=unit_nr)
313
314 cpassert(xblocks_read == natoms) ! ensure we read all blocks
315
316 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 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 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
334 TYPE(gto_basis_set_type), POINTER :: basis_set
335 TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
336 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 qs_kind_set=qs_kind_set)
341
342 IF (ikind > SIZE(atomic_kind_set) .OR. ikind > SIZE(qs_kind_set)) &
343 cpabort("Some kinds are missing.")
344
345 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, &
350 CALL pao_param_count(pao, qs_env, ikind=ikind, nparams=nparams)
351
352 IF (pao_kind%nparams /= nparams) &
353 cpabort("Number of parameters do not match")
354 IF (trim(pao_kind%name) .NE. trim(name)) &
355 cpabort("Kind names do not match")
356 IF (pao_kind%z /= z) &
357 cpabort("Atomic numbers do not match")
358 IF (trim(pao_kind%prim_basis_name) .NE. trim(basis_set%name)) &
359 cpabort("Primary Basis-set name does not match")
360 IF (pao_kind%prim_basis_size /= basis_set%nsgf) &
361 cpabort("Primary Basis-set size does not match")
362 IF (pao_kind%pao_basis_size /= pao_basis_size) &
363 cpabort("PAO basis size does not match")
364 IF (SIZE(pao_kind%pao_potentials) /= SIZE(pao_potentials)) &
365 cpabort("Number of PAO_POTENTIALS does not match")
366
367 DO ipot = 1, SIZE(pao_potentials)
368 IF (pao_kind%pao_potentials(ipot)%maxl /= pao_potentials(ipot)%maxl) &
369 cpabort("PAO_POT_MAXL does not match")
370 IF (pao_kind%pao_potentials(ipot)%max_projector /= pao_potentials(ipot)%max_projector) &
371 cpabort("PAO_POT_MAX_PROJECTOR does not match")
372 IF (pao_kind%pao_potentials(ipot)%beta /= pao_potentials(ipot)%beta) &
373 cpwarn("PAO_POT_BETA does not match")
374 IF (pao_kind%pao_potentials(ipot)%weight /= pao_potentials(ipot)%weight) &
375 cpwarn("PAO_POT_WEIGHT does not match")
376 END DO
377
378 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 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 CALL timeset(routinen, handle)
400 logger => cp_get_default_logger()
401
402 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 do_backup=.true.)
413
414 ! although just rank-0 writes the trajectory it requires collective MPI calls
415 unit_max = unit_nr
416 CALL para_env%max(unit_max)
417 IF (unit_max > 0) THEN
418 IF (pao%iw > 0) WRITE (pao%iw, '(A,A)') " PAO| Writing restart file."
419 IF (unit_nr > 0) &
420 CALL write_restart_header(pao, qs_env, energy, unit_nr)
421
422 CALL pao_write_diagonal_blocks(para_env, pao%matrix_X, "Xblock", unit_nr)
423
424 END IF
425
426 ! close file
427 IF (unit_nr > 0) WRITE (unit_nr, '(A)') "THE_END"
428 CALL cp_print_key_finished_output(unit_nr, logger, input, printkey_section)
429
430 CALL timestop(handle)
431 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 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 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
448 LOGICAL :: found
449 REAL(dp), DIMENSION(:, :), POINTER :: local_block, mpi_buffer
450
451 !TODO: this is a serial algorithm
452 CALL dbcsr_get_info(matrix, row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
453 cpassert(SIZE(row_blk_sizes) == SIZE(col_blk_sizes))
454 natoms = SIZE(row_blk_sizes)
455
456 DO iatom = 1, natoms
457 ALLOCATE (mpi_buffer(row_blk_sizes(iatom), col_blk_sizes(iatom)))
458 NULLIFY (local_block)
459 CALL dbcsr_get_block_p(matrix=matrix, row=iatom, col=iatom, block=local_block, found=found)
460 IF (ASSOCIATED(local_block)) THEN
461 IF (SIZE(local_block) > 0) & ! catch corner-case
462 mpi_buffer(:, :) = local_block(:, :)
463 ELSE
464 mpi_buffer(:, :) = 0.0_dp
465 END IF
466
467 CALL para_env%sum(mpi_buffer)
468 IF (unit_nr > 0) THEN
469 WRITE (unit_nr, fmt="(A,1X,I10,1X)", advance='no') label, iatom
470 WRITE (unit_nr, *) mpi_buffer
471 END IF
472 DEALLOCATE (mpi_buffer)
473 END DO
474
475 ! flush
476 IF (unit_nr > 0) FLUSH (unit_nr)
477
478 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 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 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 TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
500 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
501 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 qs_kind_set=qs_kind_set)
508
509 WRITE (unit_nr, "(A,5X,I0)") "Version", file_format_version
510 WRITE (unit_nr, "(A,5X,F20.10)") "Energy", energy
511 WRITE (unit_nr, "(A,5X,I0)") "Step", pao%istep
512 WRITE (unit_nr, "(A,5X,A)") "Parametrization", id2str(pao%parameterization)
513
514 ! write kinds
515 WRITE (unit_nr, "(A,5X,I0)") "Nkinds", SIZE(atomic_kind_set)
516 DO ikind = 1, SIZE(atomic_kind_set)
517 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, &
521 basis_set=basis_set)
522 CALL pao_param_count(pao, qs_env, ikind, nparams)
523 WRITE (unit_nr, "(A,5X,I10,1X,A,1X,I3)") "Kind", ikind, trim(kindname), z
524 WRITE (unit_nr, "(A,5X,I10,1X,I3)") "NParams", ikind, nparams
525 WRITE (unit_nr, "(A,5X,I10,1X,I10,1X,A)") "PrimBasis", ikind, basis_set%nsgf, trim(basis_set%name)
526 WRITE (unit_nr, "(A,5X,I10,1X,I3)") "PaoBasis", ikind, pao_basis_size
527 WRITE (unit_nr, "(A,5X,I10,1X,I3)") "NPaoPotentials", ikind, SIZE(pao_potentials)
528 DO ipot = 1, SIZE(pao_potentials)
529 WRITE (unit_nr, "(A,5X,I10,1X,I3)", advance='no') "PaoPotential", ikind, ipot
530 WRITE (unit_nr, "(1X,I3)", advance='no') pao_potentials(ipot)%maxl
531 WRITE (unit_nr, "(1X,I3)", advance='no') pao_potentials(ipot)%max_projector
532 WRITE (unit_nr, "(1X,F20.16)", advance='no') pao_potentials(ipot)%beta
533 WRITE (unit_nr, "(1X,F20.16)") pao_potentials(ipot)%weight
534 END DO
535 END DO
536
537 ! write cell
538 WRITE (unit_nr, fmt="(A,5X)", advance='no') "Cell"
539 WRITE (unit_nr, *) cell%hmat*angstrom
540
541 ! write atoms
542 WRITE (unit_nr, "(A,5X,I0)") "Natoms", SIZE(particle_set)
543 DO iatom = 1, SIZE(particle_set)
544 kindname = particle_set(iatom)%atomic_kind%name
545 WRITE (unit_nr, fmt="(A,5X,I10,5X,A,1X)", advance='no') "Atom ", iatom, trim(kindname)
546 WRITE (unit_nr, *) particle_set(iatom)%r*angstrom
547 END DO
548
549 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 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 CALL timeset(routinen, handle)
573
574 NULLIFY (dft_section)
575
576 logger => cp_get_default_logger()
577 output_unit = cp_logger_get_default_io_unit(logger)
578
579 CALL get_qs_env(qs_env, input=input)
580 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 "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 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
586
587 IF (do_ks_csr_write .AND. (.NOT. do_kpoints)) THEN
588 CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%THRESHOLD", r_val=thld)
589 CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr)
590 CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%BINARY", l_val=bin)
591
592 IF (bin) THEN
593 fileformat = "UNFORMATTED"
594 ELSE
595 fileformat = "FORMATTED"
596 END IF
597
598 DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
599
600 IF (dbcsr_has_symmetry(ls_scf_env%matrix_ks(ispin))) THEN
601 CALL dbcsr_desymmetrize(ls_scf_env%matrix_ks(ispin), matrix_ks_nosym)
602 ELSE
603 CALL dbcsr_copy(matrix_ks_nosym, ls_scf_env%matrix_ks(ispin))
604 END IF
605
606 CALL dbcsr_csr_create_from_dbcsr(matrix_ks_nosym, ks_mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
607 CALL dbcsr_convert_dbcsr_to_csr(matrix_ks_nosym, ks_mat_csr)
608
609 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 file_status="REPLACE", file_form=fileformat)
613 CALL dbcsr_csr_write(ks_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
614
615 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%KS_CSR_WRITE")
616
617 CALL dbcsr_csr_destroy(ks_mat_csr)
618 CALL dbcsr_release(matrix_ks_nosym)
619 END DO
620 END IF
621
622 CALL timestop(handle)
623
624 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 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 CALL timeset(routinen, handle)
648
649 NULLIFY (dft_section)
650
651 logger => cp_get_default_logger()
652 output_unit = cp_logger_get_default_io_unit(logger)
653
654 CALL get_qs_env(qs_env, input=input)
655 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 "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 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
661
662 IF (do_s_csr_write .AND. (.NOT. do_kpoints)) THEN
663 CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%THRESHOLD", r_val=thld)
664 CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr)
665 CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%BINARY", l_val=bin)
666
667 IF (bin) THEN
668 fileformat = "UNFORMATTED"
669 ELSE
670 fileformat = "FORMATTED"
671 END IF
672
673 IF (dbcsr_has_symmetry(ls_scf_env%matrix_s)) THEN
674 CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_s_nosym)
675 ELSE
676 CALL dbcsr_copy(matrix_s_nosym, ls_scf_env%matrix_s)
677 END IF
678
679 CALL dbcsr_csr_create_from_dbcsr(matrix_s_nosym, s_mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
680 CALL dbcsr_convert_dbcsr_to_csr(matrix_s_nosym, s_mat_csr)
681
682 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 file_status="REPLACE", file_form=fileformat)
686 CALL dbcsr_csr_write(s_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
687
688 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%S_CSR_WRITE")
689
690 CALL dbcsr_csr_destroy(s_mat_csr)
691 CALL dbcsr_release(matrix_s_nosym)
692 END IF
693
694 CALL timestop(handle)
695
696 END SUBROUTINE pao_write_s_matrix_csr
697
698END MODULE pao_io
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition cell_types.F:15
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Types needed for a linear scaling quickstep SCF run based on the density matrix.
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
Interface to the message passing library MPI.
character(len=11) function, public id2str(id)
Helper routine.
Definition pao_input.F:216
Routines for reading and writing restart files.
Definition pao_io.F:12
subroutine, public pao_kinds_ensure_equal(pao, qs_env, ikind, pao_kind)
Ensure that the kind read from the restart is equal to the kind curretly in use.
Definition pao_io.F:326
subroutine, public pao_write_restart(pao, qs_env, energy)
Writes restart file.
Definition pao_io.F:387
subroutine, public pao_write_ks_matrix_csr(qs_env, ls_scf_env)
writing the KS matrix (in terms of the PAO basis) in csr format into a file
Definition pao_io.F:558
subroutine, public pao_read_raw(filename, param, hmat, kinds, atom2kind, positions, xblocks, ml_range)
Reads a restart file into temporary datastructures.
Definition pao_io.F:183
subroutine, public pao_write_s_matrix_csr(qs_env, ls_scf_env)
writing the overlap matrix (in terms of the PAO basis) in csr format into a file
Definition pao_io.F:633
subroutine, public pao_read_restart(pao, qs_env)
Reads restart file.
Definition pao_io.F:87
Front-End for any PAO parametrization.
Definition pao_param.F:12
subroutine, public pao_param_count(pao, qs_env, ikind, nparams)
Returns the number of parameters for given atomic kind.
Definition pao_param.F:173
Factory routines for potentials used e.g. by pao_param_exp and pao_ml.
Types used by the PAO machinery.
Definition pao_types.F:12
Define the data structure for the particle information.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Holds information about a PAO potential.
Provides all information about a quickstep kind.