Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to read the binary restart file of CP2K
10 : !> \author Matthias Krack (MK)
11 : !> \par History
12 : !> - Creation (17.02.2011,MK)
13 : !> \version 1.0
14 : ! **************************************************************************************************
15 : MODULE input_cp2k_binary_restarts
16 :
17 : USE cell_types, ONLY: cell_transform_input_cartesian,&
18 : cell_type
19 : USE cp_files, ONLY: open_file
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_get_default_io_unit,&
22 : cp_logger_type,&
23 : cp_to_string
24 : USE cp_output_handling, ONLY: cp_print_key_unit_nr,&
25 : debug_print_level
26 : USE extended_system_types, ONLY: lnhc_parameters_type
27 : USE input_section_types, ONLY: section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: default_path_length,&
30 : default_string_length,&
31 : dp
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE particle_types, ONLY: particle_type
34 : USE physcon, ONLY: angstrom
35 : USE print_messages, ONLY: print_message
36 : USE string_table, ONLY: id2str,&
37 : s2s,&
38 : str2id
39 : USE topology_types, ONLY: atom_info_type,&
40 : topology_parameters_type
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_binary_restarts'
48 :
49 : PUBLIC :: read_binary_coordinates, &
50 : read_binary_cs_coordinates, &
51 : read_binary_thermostats_nose, &
52 : read_binary_velocities
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Read the input section &COORD from an external file written in
58 : !> binary format.
59 : !> \param topology ...
60 : !> \param root_section ...
61 : !> \param para_env ...
62 : !> \param subsys_section ...
63 : !> \param binary_file_read ...
64 : !> \par History
65 : !> - Creation (10.02.2011,MK)
66 : !> \author Matthias Krack (MK)
67 : !> \version 1.0
68 : ! **************************************************************************************************
69 11587 : SUBROUTINE read_binary_coordinates(topology, root_section, para_env, &
70 : subsys_section, binary_file_read)
71 :
72 : TYPE(topology_parameters_type) :: topology
73 : TYPE(section_vals_type), POINTER :: root_section
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 : TYPE(section_vals_type), POINTER :: subsys_section
76 : LOGICAL, INTENT(OUT) :: binary_file_read
77 :
78 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_binary_coordinates'
79 :
80 : CHARACTER(LEN=default_path_length) :: binary_restart_file_name
81 : CHARACTER(LEN=default_string_length) :: string
82 : INTEGER :: binary_restart_file_unit, handle, iatom, &
83 : ikind, istat, iw, natom, natomkind, &
84 : ncore, nmolecule, nmoleculekind, nshell
85 11587 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf, id_name
86 : TYPE(atom_info_type), POINTER :: atom_info
87 : TYPE(cp_logger_type), POINTER :: logger
88 :
89 11587 : CALL timeset(routineN, handle)
90 :
91 11587 : NULLIFY (logger)
92 11587 : CPASSERT(ASSOCIATED(root_section))
93 11587 : CPASSERT(ASSOCIATED(para_env))
94 11587 : CPASSERT(ASSOCIATED(subsys_section))
95 11587 : logger => cp_get_default_logger()
96 :
97 11587 : binary_file_read = .FALSE.
98 :
99 : CALL section_vals_val_get(root_section, "EXT_RESTART%BINARY_RESTART_FILE_NAME", &
100 11587 : c_val=binary_restart_file_name)
101 :
102 11587 : IF (TRIM(ADJUSTL(binary_restart_file_name)) == "") THEN
103 11541 : CALL timestop(handle)
104 11541 : RETURN
105 : END IF
106 :
107 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XYZ_INFO", &
108 46 : extension=".subsysLog")
109 :
110 46 : natomkind = 0
111 46 : natom = 0
112 46 : ncore = 0
113 46 : nshell = 0
114 46 : nmoleculekind = 0
115 46 : nmolecule = 0
116 :
117 : ! Open binary restart file and read number atomic kinds, atoms, etc.
118 46 : IF (para_env%is_source()) THEN
119 : CALL open_file(file_name=binary_restart_file_name, &
120 : file_status="OLD", &
121 : file_form="UNFORMATTED", &
122 : file_action="READWRITE", &
123 : file_position="REWIND", &
124 : unit_number=binary_restart_file_unit, &
125 23 : debug=(iw > 0))
126 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) &
127 23 : natomkind, natom, ncore, nshell, nmoleculekind, nmolecule
128 23 : IF (istat /= 0) THEN
129 : CALL stop_read("natomkind,natom,ncore,nshell,nmoleculekind,nmolecule "// &
130 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
131 0 : binary_restart_file_unit)
132 : END IF
133 23 : IF (iw > 0) THEN
134 : WRITE (UNIT=iw, FMT="(T2,A,T71,I10)") &
135 23 : "Number of atomic kinds:", natomkind, &
136 23 : "Number of atoms:", natom, &
137 23 : "Number of cores (only core-shell model):", ncore, &
138 23 : "Number of shells (only core-shell model):", nshell, &
139 23 : "Number of molecule kinds:", nmoleculekind, &
140 46 : "Number of molecules", nmolecule
141 : END IF
142 : END IF
143 :
144 46 : CALL para_env%bcast(natomkind)
145 46 : CALL para_env%bcast(natom)
146 46 : CALL para_env%bcast(ncore)
147 46 : CALL para_env%bcast(nshell)
148 46 : CALL para_env%bcast(nmoleculekind)
149 46 : CALL para_env%bcast(nmolecule)
150 :
151 138 : ALLOCATE (id_name(natomkind))
152 : ! Read atomic kind names
153 138 : DO ikind = 1, natomkind
154 92 : IF (para_env%is_source()) THEN
155 46 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) string
156 46 : IF (istat /= 0) CALL stop_read("string (IOSTAT = "// &
157 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
158 0 : binary_restart_file_unit)
159 : END IF
160 92 : CALL para_env%bcast(string)
161 138 : id_name(ikind) = str2id(string)
162 : END DO
163 :
164 : ! Allocate and initialise atom_info array
165 46 : atom_info => topology%atom_info
166 138 : ALLOCATE (atom_info%id_molname(natom))
167 4462 : atom_info%id_molname(:) = 0
168 92 : ALLOCATE (atom_info%id_resname(natom))
169 4462 : atom_info%id_resname(:) = 0
170 92 : ALLOCATE (atom_info%resid(natom))
171 4462 : atom_info%resid = 1
172 92 : ALLOCATE (atom_info%id_atmname(natom))
173 4462 : atom_info%id_atmname = 0
174 138 : ALLOCATE (atom_info%r(3, natom))
175 17710 : atom_info%r(:, :) = 0.0_dp
176 138 : ALLOCATE (atom_info%atm_mass(natom))
177 4462 : atom_info%atm_mass(:) = HUGE(0.0_dp)
178 92 : ALLOCATE (atom_info%atm_charge(natom))
179 4462 : atom_info%atm_charge(:) = -HUGE(0.0_dp)
180 92 : ALLOCATE (atom_info%occup(natom))
181 4462 : atom_info%occup(:) = 0.0_dp
182 92 : ALLOCATE (atom_info%beta(natom))
183 4462 : atom_info%beta(:) = 0.0_dp
184 92 : ALLOCATE (atom_info%id_element(natom))
185 4462 : atom_info%id_element(:) = 0
186 92 : ALLOCATE (ibuf(natom))
187 :
188 : ! Read atomic kind number of each atom
189 46 : IF (para_env%is_source()) THEN
190 23 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) ibuf(1:natom)
191 23 : IF (istat /= 0) CALL stop_read("ibuf (IOSTAT = "// &
192 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
193 0 : binary_restart_file_unit)
194 : END IF
195 46 : CALL para_env%bcast(ibuf)
196 4462 : DO iatom = 1, natom
197 4416 : ikind = ibuf(iatom)
198 4416 : atom_info%id_atmname(iatom) = id_name(ikind)
199 4462 : atom_info%id_element(iatom) = id_name(ikind)
200 : END DO
201 46 : DEALLOCATE (id_name)
202 :
203 : ! Read atomic coordinates
204 46 : IF (para_env%is_source()) THEN
205 8855 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) atom_info%r(1:3, 1:natom)
206 23 : IF (istat /= 0) CALL stop_read("atom_info%r(1:3,1:natom) (IOSTAT = "// &
207 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
208 0 : binary_restart_file_unit)
209 : END IF
210 35374 : CALL para_env%bcast(atom_info%r)
211 46 : IF (ASSOCIATED(topology%cell_muc)) THEN
212 4462 : DO iatom = 1, natom
213 4462 : CALL cell_transform_input_cartesian(topology%cell_muc, atom_info%r(:, iatom))
214 : END DO
215 : END IF
216 :
217 : ! Read molecule information if available
218 46 : IF (nmolecule > 0) THEN
219 138 : ALLOCATE (id_name(nmoleculekind))
220 : ! Read molecule kind names
221 92 : DO ikind = 1, nmoleculekind
222 46 : IF (para_env%is_source()) THEN
223 23 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) string
224 23 : IF (istat /= 0) CALL stop_read("string (IOSTAT = "// &
225 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
226 0 : binary_restart_file_unit)
227 : END IF
228 46 : CALL para_env%bcast(string)
229 92 : id_name(ikind) = str2id(string)
230 : END DO
231 : ! Read molecule kind numbers
232 46 : IF (para_env%is_source()) THEN
233 23 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) ibuf(1:natom)
234 23 : IF (istat /= 0) CALL stop_read("ibuf(1:natom) (IOSTAT = "// &
235 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
236 0 : binary_restart_file_unit)
237 : END IF
238 46 : CALL para_env%bcast(ibuf)
239 4462 : DO iatom = 1, natom
240 4416 : ikind = ibuf(iatom)
241 4462 : atom_info%id_molname(iatom) = id_name(ikind)
242 : END DO
243 46 : DEALLOCATE (id_name)
244 : ! Read molecule index which is used also as residue id
245 46 : IF (para_env%is_source()) THEN
246 2231 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) atom_info%resid(1:natom)
247 23 : IF (istat /= 0) CALL stop_read("atom_info%resid(1:natom) (IOSTAT = "// &
248 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
249 0 : binary_restart_file_unit)
250 : END IF
251 8878 : CALL para_env%bcast(atom_info%resid)
252 4462 : DO iatom = 1, natom
253 4462 : atom_info%id_resname(iatom) = str2id(s2s(cp_to_string(atom_info%resid(iatom))))
254 : END DO
255 : END IF
256 46 : DEALLOCATE (ibuf)
257 :
258 : !MK to be checked ...
259 46 : topology%aa_element = .TRUE.
260 46 : topology%molname_generated = .FALSE.
261 46 : topology%natoms = natom
262 :
263 46 : IF (iw > 0) THEN
264 : WRITE (UNIT=iw, FMT="(T2,A)") &
265 : "BEGIN of COORD section data [Angstrom] read in binary format from file "// &
266 23 : TRIM(binary_restart_file_name)
267 2231 : DO iatom = 1, natom
268 : WRITE (UNIT=iw, FMT="(T2,A2,3(1X,ES25.16),2(1X,A))") &
269 2208 : TRIM(ADJUSTL(id2str(atom_info%id_atmname(iatom)))), &
270 8832 : atom_info%r(1:3, iatom)*angstrom, &
271 2208 : TRIM(ADJUSTL(id2str(atom_info%id_molname(iatom)))), &
272 4439 : TRIM(ADJUSTL(id2str(atom_info%id_resname(iatom))))
273 : END DO
274 : WRITE (UNIT=iw, FMT="(T2,A)") &
275 : "END of COORD section data [Angstrom] read from binary restart file "// &
276 23 : TRIM(binary_restart_file_name)
277 : END IF
278 :
279 : ! Keep the file connected: the next reader continues at the current position
280 :
281 46 : binary_file_read = .TRUE.
282 :
283 46 : CALL timestop(handle)
284 :
285 11587 : END SUBROUTINE read_binary_coordinates
286 :
287 : ! **************************************************************************************************
288 : !> \brief Read the input section &CORE_COORD or &SHELL_COORD from an external
289 : !> file written in binary format.
290 : !> \param prefix ...
291 : !> \param particle_set ...
292 : !> \param root_section ...
293 : !> \param subsys_section ...
294 : !> \param binary_file_read ...
295 : !> \param cell ...
296 : !> \par History
297 : !> - Creation (17.02.2011,MK)
298 : !> \author Matthias Krack (MK)
299 : !> \version 1.0
300 : ! **************************************************************************************************
301 5398 : SUBROUTINE read_binary_cs_coordinates(prefix, particle_set, root_section, &
302 : subsys_section, binary_file_read, cell)
303 :
304 : CHARACTER(LEN=*), INTENT(IN) :: prefix
305 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
306 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
307 : LOGICAL, INTENT(OUT) :: binary_file_read
308 : TYPE(cell_type), OPTIONAL, POINTER :: cell
309 :
310 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_binary_cs_coordinates'
311 :
312 : CHARACTER(LEN=default_path_length) :: binary_restart_file_name, message
313 : CHARACTER(LEN=default_string_length) :: section_label, section_name
314 : INTEGER :: binary_restart_file_unit, handle, &
315 : iparticle, istat, iw, nbuf, nparticle
316 5398 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf
317 : LOGICAL :: exit_routine
318 5398 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rbuf
319 : TYPE(cp_logger_type), POINTER :: logger
320 : TYPE(mp_para_env_type), POINTER :: para_env
321 :
322 5398 : CALL timeset(routineN, handle)
323 :
324 5398 : NULLIFY (logger)
325 5398 : CPASSERT(ASSOCIATED(root_section))
326 5398 : CPASSERT(ASSOCIATED(subsys_section))
327 5398 : logger => cp_get_default_logger()
328 5398 : para_env => logger%para_env
329 :
330 5398 : binary_file_read = .FALSE.
331 :
332 5398 : IF (ASSOCIATED(particle_set)) THEN
333 504 : exit_routine = .FALSE.
334 504 : nparticle = SIZE(particle_set)
335 : ELSE
336 4894 : exit_routine = .TRUE.
337 4894 : nparticle = 0
338 : END IF
339 :
340 : CALL section_vals_val_get(root_section, "EXT_RESTART%BINARY_RESTART_FILE_NAME", &
341 5398 : c_val=binary_restart_file_name)
342 :
343 5398 : IF (TRIM(ADJUSTL(binary_restart_file_name)) == "") THEN
344 5306 : CALL timestop(handle)
345 5306 : RETURN
346 : END IF
347 :
348 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XYZ_INFO", &
349 92 : extension=".subsysLog")
350 :
351 92 : section_name = prefix//" COORDINATES"
352 :
353 : ! Open binary restart file at last position
354 92 : IF (para_env%is_source()) THEN
355 : CALL open_file(file_name=TRIM(binary_restart_file_name), &
356 : file_status="OLD", &
357 : file_form="UNFORMATTED", &
358 : file_action="READWRITE", &
359 : file_position="ASIS", &
360 : unit_number=binary_restart_file_unit, &
361 46 : debug=(iw > 0))
362 46 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) section_label, nbuf
363 46 : IF (istat /= 0) CALL stop_read("section_label, nbuf -> "//TRIM(section_label)//", "// &
364 : TRIM(ADJUSTL(cp_to_string(nbuf)))// &
365 : " (IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//"). "// &
366 : "Section "//TRIM(ADJUSTL(section_name))//" was expected.", &
367 0 : binary_restart_file_unit)
368 46 : IF (TRIM(section_label) == TRIM(section_name)) THEN
369 46 : IF (nbuf /= nparticle) THEN
370 16 : IF (iw > 0) THEN
371 : message = "INFO: The requested number of "//TRIM(section_name)//" ("// &
372 : TRIM(ADJUSTL(cp_to_string(nparticle)))//") does not match the "// &
373 : "number ("//TRIM(ADJUSTL(cp_to_string(nbuf)))//") available from the "// &
374 : "binary restart file <"//TRIM(binary_restart_file_name)// &
375 16 : ">. The restart file information is ignored."
376 16 : CALL print_message(message, iw, 1, 1, 1)
377 : END IF
378 : ! Ignore this section
379 16 : IF (nbuf > 0) THEN
380 : ! Perform dummy read
381 24 : ALLOCATE (rbuf(3, nbuf))
382 8 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:3, 1:nbuf)
383 8 : IF (istat /= 0) CALL stop_read("rbuf(1:3,1:nbuf) -> "//prefix// &
384 : " coordinates (IOSTAT = "// &
385 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
386 0 : binary_restart_file_unit)
387 8 : DEALLOCATE (rbuf)
388 24 : ALLOCATE (ibuf(nbuf))
389 8 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) ibuf(1:nbuf)
390 8 : IF (istat /= 0) CALL stop_read("ibuf(1:nparticle) -> atomic indices of the "// &
391 : TRIM(section_name)//" (IOSTAT = "// &
392 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
393 0 : binary_restart_file_unit)
394 8 : DEALLOCATE (ibuf)
395 : END IF
396 16 : exit_routine = .TRUE.
397 : ELSE
398 30 : IF (iw > 0) THEN
399 : WRITE (UNIT=iw, FMT="(T2,A,T71,I10)") &
400 30 : "Number of "//prefix//" particles:", nparticle
401 : END IF
402 30 : IF (nparticle == 0) exit_routine = .TRUE.
403 : END IF
404 : ELSE
405 : CALL cp_abort(__LOCATION__, &
406 : "Section label <"//TRIM(section_label)//"> read from the "// &
407 : "binary restart file <"//TRIM(binary_restart_file_name)// &
408 : "> does not match the requested section name <"// &
409 0 : TRIM(section_name)//">.")
410 : END IF
411 : END IF
412 :
413 92 : CALL para_env%bcast(exit_routine)
414 92 : IF (exit_routine) THEN
415 : ! Keep the file connected: the next reader continues at the current position
416 60 : CALL timestop(handle)
417 60 : RETURN
418 : END IF
419 :
420 32 : CPASSERT(nparticle > 0)
421 :
422 96 : ALLOCATE (rbuf(3, nparticle))
423 :
424 32 : IF (para_env%is_source()) THEN
425 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:3, 1:nparticle)
426 16 : IF (istat /= 0) CALL stop_read("rbuf(1:3,1:nparticle) -> "//prefix// &
427 : " coordinates (IOSTAT = "// &
428 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
429 0 : binary_restart_file_unit)
430 : END IF
431 32 : CALL para_env%bcast(rbuf)
432 :
433 3104 : DO iparticle = 1, nparticle
434 12288 : particle_set(iparticle)%r(1:3) = rbuf(1:3, iparticle)
435 3104 : IF (PRESENT(cell)) THEN
436 3072 : IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, particle_set(iparticle)%r(1:3))
437 : END IF
438 : END DO
439 :
440 32 : DEALLOCATE (rbuf)
441 :
442 96 : ALLOCATE (ibuf(nparticle))
443 :
444 32 : IF (para_env%is_source()) THEN
445 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) ibuf(1:nparticle)
446 16 : IF (istat /= 0) CALL stop_read("ibuf(1:nparticle) -> atomic indices of the "// &
447 : TRIM(section_name)//" (IOSTAT = "// &
448 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
449 0 : binary_restart_file_unit)
450 : END IF
451 :
452 32 : CALL para_env%bcast(ibuf)
453 :
454 3104 : DO iparticle = 1, nparticle
455 3104 : particle_set(iparticle)%atom_index = ibuf(iparticle)
456 : END DO
457 :
458 32 : DEALLOCATE (ibuf)
459 :
460 32 : IF (iw > 0) THEN
461 : WRITE (UNIT=iw, FMT="(T2,A)") &
462 : "BEGIN of "//TRIM(ADJUSTL(section_name))// &
463 : " section data [Angstrom] read in binary format from file "// &
464 16 : TRIM(binary_restart_file_name)
465 1552 : DO iparticle = 1, nparticle
466 : WRITE (UNIT=iw, FMT="(T2,A2,3(1X,ES25.16),1X,I0)") &
467 1536 : TRIM(ADJUSTL(particle_set(iparticle)%atomic_kind%name)), &
468 6144 : particle_set(iparticle)%r(1:3)*angstrom, &
469 3088 : particle_set(iparticle)%atom_index
470 : END DO
471 : WRITE (UNIT=iw, FMT="(T2,A)") &
472 : "END of "//TRIM(ADJUSTL(section_name))// &
473 : " section data [Angstrom] read from binary restart file "// &
474 16 : TRIM(binary_restart_file_name)
475 : END IF
476 :
477 : ! Keep the file connected: the next reader continues at the current position
478 :
479 32 : binary_file_read = .TRUE.
480 :
481 32 : CALL timestop(handle)
482 :
483 10796 : END SUBROUTINE read_binary_cs_coordinates
484 :
485 : ! **************************************************************************************************
486 : !> \brief Read the input section &VELOCITY, &CORE_VELOCITY, or
487 : !> &SHELL_VELOCITY from an external file written in binary format.
488 : !> \param prefix ...
489 : !> \param particle_set ...
490 : !> \param root_section ...
491 : !> \param para_env ...
492 : !> \param subsys_section ...
493 : !> \param binary_file_read ...
494 : !> \param cell ...
495 : !> \par History
496 : !> - Creation (17.02.2011,MK)
497 : !> \author Matthias Krack (MK)
498 : !> \version 1.0
499 : ! **************************************************************************************************
500 5355 : SUBROUTINE read_binary_velocities(prefix, particle_set, root_section, para_env, &
501 : subsys_section, binary_file_read, cell)
502 :
503 : CHARACTER(LEN=*), INTENT(IN) :: prefix
504 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
505 : TYPE(section_vals_type), POINTER :: root_section
506 : TYPE(mp_para_env_type), POINTER :: para_env
507 : TYPE(section_vals_type), POINTER :: subsys_section
508 : LOGICAL, INTENT(OUT) :: binary_file_read
509 : TYPE(cell_type), OPTIONAL, POINTER :: cell
510 :
511 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_binary_velocities'
512 :
513 : CHARACTER(LEN=default_path_length) :: binary_restart_file_name, message
514 : CHARACTER(LEN=default_string_length) :: section_label, section_name
515 : INTEGER :: binary_restart_file_unit, handle, i, &
516 : iparticle, istat, iw, nbuf, nparticle
517 : LOGICAL :: have_velocities
518 5355 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rbuf
519 : TYPE(cp_logger_type), POINTER :: logger
520 :
521 5355 : CALL timeset(routineN, handle)
522 :
523 5355 : NULLIFY (logger)
524 5355 : CPASSERT(ASSOCIATED(root_section))
525 5355 : CPASSERT(ASSOCIATED(para_env))
526 5355 : CPASSERT(ASSOCIATED(subsys_section))
527 5355 : logger => cp_get_default_logger()
528 :
529 5355 : binary_file_read = .FALSE.
530 :
531 : CALL section_vals_val_get(root_section, "EXT_RESTART%BINARY_RESTART_FILE_NAME", &
532 5355 : c_val=binary_restart_file_name)
533 :
534 5355 : IF (TRIM(ADJUSTL(binary_restart_file_name)) == "") THEN
535 5217 : CALL timestop(handle)
536 5217 : RETURN
537 : END IF
538 :
539 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XYZ_INFO", &
540 138 : extension=".subsysLog")
541 :
542 138 : IF (LEN_TRIM(prefix) == 0) THEN
543 46 : section_name = "VELOCITIES"
544 : ELSE
545 92 : section_name = prefix//" VELOCITIES"
546 : END IF
547 :
548 138 : have_velocities = .FALSE.
549 :
550 138 : IF (ASSOCIATED(particle_set)) THEN
551 94 : nparticle = SIZE(particle_set)
552 : ELSE
553 44 : nparticle = 0
554 : END IF
555 :
556 : ! Open binary restart file at last position and check if there are
557 : ! velocities available
558 138 : IF (para_env%is_source()) THEN
559 : CALL open_file(file_name=binary_restart_file_name, &
560 : file_status="OLD", &
561 : file_form="UNFORMATTED", &
562 : file_action="READWRITE", &
563 : file_position="ASIS", &
564 : unit_number=binary_restart_file_unit, &
565 69 : debug=(iw > 0))
566 : DO
567 96 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) section_label, nbuf
568 96 : IF (istat /= 0) CALL stop_read("section_label, nbuf -> "//TRIM(section_label)//", "// &
569 : TRIM(ADJUSTL(cp_to_string(nbuf)))// &
570 : " (IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//"). "// &
571 : "Section "//TRIM(ADJUSTL(section_name))//" was expected.", &
572 0 : binary_restart_file_unit)
573 96 : IF (INDEX(section_label, "THERMOSTAT") > 0) THEN
574 27 : IF (nbuf > 0) THEN
575 : ! Ignore thermostat information
576 12 : ALLOCATE (rbuf(nbuf, 1))
577 : ! Perform dummy read
578 20 : DO i = 1, 4
579 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:nbuf, 1)
580 16 : IF (istat /= 0) CALL stop_read("rbuf(1:nbuf,1) -> "// &
581 : TRIM(ADJUSTL(section_label))// &
582 : " (IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
583 4 : binary_restart_file_unit)
584 : END DO
585 4 : DEALLOCATE (rbuf)
586 4 : IF (iw > 0) THEN
587 : message = "INFO: Ignoring section <"//TRIM(ADJUSTL(section_label))// &
588 4 : "> from binary restart file <"//TRIM(binary_restart_file_name)//">."
589 4 : CALL print_message(message, iw, 1, 1, 1)
590 : END IF
591 : END IF
592 : CYCLE
593 69 : ELSE IF (INDEX(section_label, "VELOCIT") == 0) THEN
594 : CALL cp_abort(__LOCATION__, &
595 : "Section label <"//TRIM(section_label)//"> read from the "// &
596 : "binary restart file <"//TRIM(binary_restart_file_name)// &
597 : "> does not match the requested section name <"// &
598 0 : TRIM(section_name)//">.")
599 : ELSE
600 69 : IF (nbuf > 0) have_velocities = .TRUE.
601 : EXIT
602 : END IF
603 : END DO
604 : END IF
605 :
606 138 : CALL para_env%bcast(nbuf)
607 138 : CALL para_env%bcast(have_velocities)
608 :
609 138 : IF (have_velocities) THEN
610 :
611 282 : ALLOCATE (rbuf(3, nbuf))
612 :
613 94 : IF (para_env%is_source()) THEN
614 47 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:3, 1:nbuf)
615 47 : IF (istat /= 0) CALL stop_read("rbuf(1:3,1:nbuf) -> "// &
616 : TRIM(ADJUSTL(section_name))// &
617 : " (IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
618 0 : binary_restart_file_unit)
619 : END IF
620 :
621 94 : IF (nbuf == nparticle) THEN
622 78 : CALL para_env%bcast(rbuf)
623 7566 : DO iparticle = 1, nparticle
624 29952 : particle_set(iparticle)%v(1:3) = rbuf(1:3, iparticle)
625 7566 : IF (PRESENT(cell)) THEN
626 7488 : IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, particle_set(iparticle)%v(1:3))
627 : END IF
628 : END DO
629 : ELSE
630 16 : IF (iw > 0) THEN
631 : message = "INFO: The requested number of "//TRIM(ADJUSTL(section_name))// &
632 : " ("//TRIM(ADJUSTL(cp_to_string(nparticle)))//") does not match the "// &
633 : "number ("//TRIM(ADJUSTL(cp_to_string(nbuf)))//") available from the "// &
634 : "binary restart file <"//TRIM(binary_restart_file_name)// &
635 8 : ">. The restart file information is ignored."
636 8 : CALL print_message(message, iw, 1, 1, 1)
637 : END IF
638 : END IF
639 :
640 94 : DEALLOCATE (rbuf)
641 :
642 : END IF
643 :
644 138 : IF (nbuf == nparticle) THEN
645 106 : IF (iw > 0) THEN
646 : WRITE (UNIT=iw, FMT="(T2,A)") &
647 : "BEGIN of "//TRIM(ADJUSTL(section_name))// &
648 : " section data [a.u.] read in binary format from file "// &
649 53 : TRIM(binary_restart_file_name)
650 53 : IF (have_velocities) THEN
651 3783 : DO iparticle = 1, nparticle
652 : WRITE (UNIT=iw, FMT="(T2,A2,3(1X,ES25.16))") &
653 3744 : TRIM(ADJUSTL(particle_set(iparticle)%atomic_kind%name)), &
654 18759 : particle_set(iparticle)%v(1:3)
655 : END DO
656 : ELSE
657 : WRITE (UNIT=iw, FMT="(A)") &
658 14 : "# No "//TRIM(ADJUSTL(section_name))//" available"
659 : END IF
660 : WRITE (UNIT=iw, FMT="(T2,A)") &
661 : "END of "//TRIM(ADJUSTL(section_name))// &
662 : " section data [a.u.] read from binary restart file "// &
663 53 : TRIM(binary_restart_file_name)
664 : END IF
665 106 : binary_file_read = .TRUE.
666 : END IF
667 :
668 : ! Keep the file connected: the next reader continues at the current position
669 :
670 138 : CALL timestop(handle)
671 :
672 138 : END SUBROUTINE read_binary_velocities
673 :
674 : ! **************************************************************************************************
675 : !> \brief Read the input section &THERMOSTAT for Nose thermostats from an
676 : !> external file written in binary format.
677 : !> \param prefix ...
678 : !> \param nhc ...
679 : !> \param binary_restart_file_name ...
680 : !> \param restart ...
681 : !> \param para_env ...
682 : !> \par History
683 : !> - Creation (28.02.2011,MK)
684 : !> \author Matthias Krack (MK)
685 : !> \version 1.0
686 : ! **************************************************************************************************
687 38 : SUBROUTINE read_binary_thermostats_nose(prefix, nhc, binary_restart_file_name, &
688 : restart, para_env)
689 :
690 : CHARACTER(LEN=*), INTENT(IN) :: prefix
691 : TYPE(lnhc_parameters_type), POINTER :: nhc
692 : CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name
693 : LOGICAL, INTENT(OUT) :: restart
694 : TYPE(mp_para_env_type), POINTER :: para_env
695 :
696 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_binary_thermostats_nose'
697 :
698 : CHARACTER(LEN=default_string_length) :: section_label, section_name
699 : INTEGER :: binary_restart_file_unit, handle, i, &
700 : idx, istat, j, nhc_size, output_unit
701 : LOGICAL :: debug
702 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rbuf
703 : TYPE(cp_logger_type), POINTER :: logger
704 :
705 38 : CALL timeset(routineN, handle)
706 :
707 38 : CPASSERT(ASSOCIATED(nhc))
708 38 : CPASSERT(ASSOCIATED(para_env))
709 :
710 : ! Set to .TRUE. for debug mode, i.e. all data read are written to stdout
711 38 : NULLIFY (logger)
712 38 : logger => cp_get_default_logger()
713 38 : output_unit = cp_logger_get_default_io_unit(logger)
714 :
715 38 : IF (logger%iter_info%print_level >= debug_print_level) THEN
716 12 : debug = .TRUE.
717 : ELSE
718 26 : debug = .FALSE.
719 : END IF
720 :
721 38 : restart = .FALSE.
722 :
723 38 : section_name = prefix//" THERMOSTATS"
724 :
725 : ! Open binary restart file at last position
726 38 : IF (para_env%is_source()) THEN
727 : CALL open_file(file_name=binary_restart_file_name, &
728 : file_status="OLD", &
729 : file_form="UNFORMATTED", &
730 : file_action="READWRITE", &
731 : file_position="ASIS", &
732 : unit_number=binary_restart_file_unit, &
733 19 : debug=debug)
734 19 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) section_label, nhc_size
735 19 : IF (istat /= 0) CALL stop_read("nhc_size (IOSTAT = "// &
736 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
737 0 : binary_restart_file_unit)
738 19 : IF (INDEX(section_label, "THERMOSTAT") == 0) THEN
739 : CALL cp_abort(__LOCATION__, &
740 : "Section label <"//TRIM(section_label)//"> read from the "// &
741 : "binary restart file <"//TRIM(binary_restart_file_name)// &
742 : "> does not match the requested section name <"// &
743 0 : TRIM(section_name)//">.")
744 : END IF
745 19 : IF (debug .AND. output_unit > 0) THEN
746 : WRITE (UNIT=output_unit, FMT="(T2,A,/,T2,A,I0)") &
747 : "BEGIN of "//TRIM(ADJUSTL(section_label))// &
748 : " section data read in binary format from file "// &
749 6 : TRIM(binary_restart_file_name), &
750 12 : "# nhc_size = ", nhc_size
751 : END IF
752 : END IF
753 :
754 38 : CALL para_env%bcast(nhc_size)
755 :
756 38 : IF (nhc_size > 0) THEN
757 :
758 96 : ALLOCATE (rbuf(nhc_size))
759 32 : rbuf(:) = 0.0_dp
760 :
761 : ! Read NHC section &COORD
762 32 : IF (para_env%is_source()) THEN
763 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:nhc_size)
764 16 : IF (istat /= 0) CALL stop_read("eta -> rbuf (IOSTAT = "// &
765 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
766 0 : binary_restart_file_unit)
767 16 : IF (debug .AND. output_unit > 0) THEN
768 : WRITE (UNIT=output_unit, FMT="(T2,A,/,(4(1X,ES25.16)))") &
769 4 : "&COORD", rbuf(1:nhc_size)
770 : END IF
771 : END IF
772 32 : CALL para_env%bcast(rbuf)
773 4640 : DO i = 1, SIZE(nhc%nvt, 2)
774 4608 : idx = (nhc%map_info%index(i) - 1)*nhc%nhc_len
775 18464 : DO j = 1, SIZE(nhc%nvt, 1)
776 13824 : idx = idx + 1
777 18432 : nhc%nvt(j, i)%eta = rbuf(idx)
778 : END DO
779 : END DO
780 :
781 : ! Read NHC section &VELOCITY
782 32 : IF (para_env%is_source()) THEN
783 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:nhc_size)
784 16 : IF (istat /= 0) CALL stop_read("veta -> rbuf (IOSTAT = "// &
785 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
786 0 : binary_restart_file_unit)
787 16 : IF (debug .AND. output_unit > 0) THEN
788 : WRITE (UNIT=output_unit, FMT="(T2,A,/,(4(1X,ES25.16)))") &
789 4 : "&VELOCITY", rbuf(1:nhc_size)
790 : END IF
791 : END IF
792 32 : CALL para_env%bcast(rbuf)
793 4640 : DO i = 1, SIZE(nhc%nvt, 2)
794 4608 : idx = (nhc%map_info%index(i) - 1)*nhc%nhc_len
795 18464 : DO j = 1, SIZE(nhc%nvt, 1)
796 13824 : idx = idx + 1
797 18432 : nhc%nvt(j, i)%v = rbuf(idx)
798 : END DO
799 : END DO
800 :
801 : ! Read NHC section &MASS
802 32 : IF (para_env%is_source()) THEN
803 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:nhc_size)
804 16 : IF (istat /= 0) CALL stop_read("mnhc -> rbuf (IOSTAT = "// &
805 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
806 0 : binary_restart_file_unit)
807 16 : IF (debug .AND. output_unit > 0) THEN
808 : WRITE (UNIT=output_unit, FMT="(T2,A,/,(4(1X,ES25.16)))") &
809 4 : "&MASS:", rbuf(1:nhc_size)
810 : END IF
811 : END IF
812 32 : CALL para_env%bcast(rbuf)
813 4640 : DO i = 1, SIZE(nhc%nvt, 2)
814 4608 : idx = (nhc%map_info%index(i) - 1)*nhc%nhc_len
815 18464 : DO j = 1, SIZE(nhc%nvt, 1)
816 13824 : idx = idx + 1
817 18432 : nhc%nvt(j, i)%mass = rbuf(idx)
818 : END DO
819 : END DO
820 :
821 : ! Read NHC section &FORCE
822 32 : IF (para_env%is_source()) THEN
823 16 : READ (UNIT=binary_restart_file_unit, IOSTAT=istat) rbuf(1:nhc_size)
824 16 : IF (istat /= 0) CALL stop_read("fnhc -> rbuf (IOSTAT = "// &
825 : TRIM(ADJUSTL(cp_to_string(istat)))//")", &
826 0 : binary_restart_file_unit)
827 16 : IF (debug .AND. output_unit > 0) THEN
828 : WRITE (UNIT=output_unit, FMT="(T2,A,/,(4(1X,ES25.16)))") &
829 4 : "&FORCE", rbuf(1:nhc_size)
830 : END IF
831 : END IF
832 32 : CALL para_env%bcast(rbuf)
833 4640 : DO i = 1, SIZE(nhc%nvt, 2)
834 4608 : idx = (nhc%map_info%index(i) - 1)*nhc%nhc_len
835 18464 : DO j = 1, SIZE(nhc%nvt, 1)
836 13824 : idx = idx + 1
837 18432 : nhc%nvt(j, i)%f = rbuf(idx)
838 : END DO
839 : END DO
840 :
841 32 : DEALLOCATE (rbuf)
842 :
843 32 : restart = .TRUE.
844 :
845 : END IF
846 :
847 38 : IF (para_env%is_source()) THEN
848 19 : IF (debug .AND. output_unit > 0) THEN
849 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
850 : "END of"//TRIM(ADJUSTL(section_label))// &
851 : " section data read in binary format from file "// &
852 6 : TRIM(binary_restart_file_name)
853 : END IF
854 : ! Keep the file connected: the next reader continues at the current position
855 : END IF
856 :
857 38 : CALL timestop(handle)
858 :
859 38 : END SUBROUTINE read_binary_thermostats_nose
860 :
861 : ! **************************************************************************************************
862 : !> \brief Print an error message and stop the program execution in case of a
863 : !> read error.
864 : !> \param object ...
865 : !> \param unit_number ...
866 : !> \par History
867 : !> - Creation (15.02.2011,MK)
868 : !> \author Matthias Krack (MK)
869 : !> \note
870 : !> object : Name of the data object for which I/O operation failed
871 : !> unit_number: Logical unit number of the file read from
872 : ! **************************************************************************************************
873 0 : SUBROUTINE stop_read(object, unit_number)
874 : CHARACTER(LEN=*), INTENT(IN) :: object
875 : INTEGER, INTENT(IN) :: unit_number
876 :
877 : CHARACTER(LEN=2*default_path_length) :: message
878 : CHARACTER(LEN=default_path_length) :: file_name
879 : LOGICAL :: file_exists
880 :
881 0 : IF (unit_number >= 0) THEN
882 0 : INQUIRE (UNIT=unit_number, EXIST=file_exists)
883 : ELSE
884 0 : file_exists = .FALSE.
885 : END IF
886 0 : IF (file_exists) THEN
887 0 : INQUIRE (UNIT=unit_number, NAME=file_name)
888 : WRITE (UNIT=message, FMT="(A)") &
889 : "An error occurred reading data object <"//TRIM(ADJUSTL(object))// &
890 0 : "> from file <"//TRIM(ADJUSTL(file_name))//">"
891 : ELSE
892 : WRITE (UNIT=message, FMT="(A,I0,A)") &
893 : "Could not read data object <"//TRIM(ADJUSTL(object))// &
894 0 : "> from logical unit ", unit_number, ". The I/O unit does not exist."
895 : END IF
896 :
897 0 : CPABORT(message)
898 :
899 0 : END SUBROUTINE stop_read
900 :
901 : END MODULE input_cp2k_binary_restarts
|
