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