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