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 Set of routines to dump the restart file of CP2K
10 : !> \par History
11 : !> 01.2006 [created] Teodoro Laino
12 : ! **************************************************************************************************
13 : MODULE input_cp2k_restarts
14 :
15 : USE al_system_types, ONLY: al_system_type
16 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
17 : USE averages_types, ONLY: average_quantities_type
18 : USE cp2k_info, ONLY: write_restart_header
19 : USE cp_linked_list_input, ONLY: cp_sll_val_create,&
20 : cp_sll_val_get_length,&
21 : cp_sll_val_type
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_get_default_io_unit,&
24 : cp_logger_type,&
25 : cp_to_string
26 : USE cp_output_handling, ONLY: cp_p_file,&
27 : cp_print_key_finished_output,&
28 : cp_print_key_should_output,&
29 : cp_print_key_unit_nr
30 : USE cp_subsys_types, ONLY: cp_subsys_get,&
31 : cp_subsys_type
32 : USE csvr_system_types, ONLY: csvr_system_type
33 : USE extended_system_types, ONLY: lnhc_parameters_type,&
34 : map_info_type,&
35 : npt_info_type
36 : USE force_env_types, ONLY: force_env_get,&
37 : force_env_type,&
38 : multiple_fe_list
39 : USE gle_system_types, ONLY: gle_type
40 : USE helium_types, ONLY: helium_solvent_p_type
41 : USE input_constants, ONLY: &
42 : do_band_collective, do_thermo_al, do_thermo_csvr, do_thermo_gle, &
43 : do_thermo_no_communication, do_thermo_nose, mol_dyn_run, mon_car_run, pint_run
44 : USE input_restart_force_eval, ONLY: update_force_eval
45 : USE input_restart_rng, ONLY: section_rng_val_set
46 : USE input_section_types, ONLY: &
47 : section_get_keyword_index, section_type, section_vals_add_values, section_vals_get, &
48 : section_vals_get_subs_vals, section_vals_get_subs_vals3, section_vals_remove_values, &
49 : section_vals_type, section_vals_val_get, section_vals_val_set, section_vals_val_unset, &
50 : section_vals_write
51 : USE input_val_types, ONLY: val_create,&
52 : val_release,&
53 : val_type
54 : USE kinds, ONLY: default_path_length,&
55 : default_string_length,&
56 : dp,&
57 : dp_size,&
58 : int_size
59 : USE md_environment_types, ONLY: get_md_env,&
60 : md_environment_type
61 : USE memory_utilities, ONLY: reallocate
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE metadynamics_types, ONLY: meta_env_type
64 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
65 : USE molecule_list_types, ONLY: molecule_list_type
66 : USE neb_types, ONLY: neb_var_type
67 : USE parallel_rng_types, ONLY: rng_record_length
68 : USE particle_list_types, ONLY: particle_list_type
69 : USE particle_types, ONLY: get_particle_pos_or_vel,&
70 : particle_type
71 : USE physcon, ONLY: angstrom
72 : USE pint_transformations, ONLY: pint_u2x
73 : USE pint_types, ONLY: pint_env_type,&
74 : thermostat_gle,&
75 : thermostat_nose,&
76 : thermostat_piglet,&
77 : thermostat_pile,&
78 : thermostat_qtb
79 : USE simpar_types, ONLY: simpar_type
80 : USE string_utilities, ONLY: string_to_ascii
81 : USE thermostat_types, ONLY: thermostat_type
82 : USE thermostat_utils, ONLY: communication_thermo_low2,&
83 : get_kin_energies
84 : #include "../base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_restarts'
91 :
92 : PUBLIC :: write_restart
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief checks if a restart needs to be written and does so, updating all necessary fields
98 : !> in the input file. This is a relatively simple wrapper routine.
99 : !> \param md_env ...
100 : !> \param force_env ...
101 : !> \param root_section ...
102 : !> \param coords ...
103 : !> \param vels ...
104 : !> \param pint_env ...
105 : !> \param helium_env ...
106 : !> \par History
107 : !> 03.2006 created [Joost VandeVondele]
108 : !> \author Joost VandeVondele
109 : ! **************************************************************************************************
110 103836 : SUBROUTINE write_restart(md_env, force_env, root_section, &
111 : coords, vels, pint_env, helium_env)
112 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
113 : TYPE(force_env_type), OPTIONAL, POINTER :: force_env
114 : TYPE(section_vals_type), POINTER :: root_section
115 : TYPE(neb_var_type), OPTIONAL, POINTER :: coords, vels
116 : TYPE(pint_env_type), INTENT(IN), OPTIONAL :: pint_env
117 : TYPE(helium_solvent_p_type), DIMENSION(:), &
118 : OPTIONAL, POINTER :: helium_env
119 :
120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_restart'
121 : CHARACTER(LEN=30), DIMENSION(2), PARAMETER :: &
122 : keys = ["PRINT%RESTART_HISTORY", "PRINT%RESTART "]
123 :
124 : INTEGER :: handle, ikey, ires, log_unit, nforce_eval
125 : LOGICAL :: save_mem, write_binary_restart_file
126 : TYPE(cp_logger_type), POINTER :: logger
127 : TYPE(section_vals_type), POINTER :: global_section, motion_section, sections
128 :
129 51918 : CALL timeset(routineN, handle)
130 :
131 51918 : logger => cp_get_default_logger()
132 51918 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
133 :
134 51918 : NULLIFY (global_section)
135 51918 : global_section => section_vals_get_subs_vals(root_section, "GLOBAL")
136 51918 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
137 :
138 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
139 51918 : motion_section, keys(1)), cp_p_file) .OR. &
140 : BTEST(cp_print_key_should_output(logger%iter_info, &
141 : motion_section, keys(2)), cp_p_file)) THEN
142 :
143 14046 : sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
144 14046 : CALL section_vals_get(sections, n_repetition=nforce_eval)
145 : CALL section_vals_val_get(motion_section, "PRINT%RESTART%SPLIT_RESTART_FILE", &
146 14046 : l_val=write_binary_restart_file)
147 :
148 14046 : IF (write_binary_restart_file) THEN
149 136 : CALL update_subsys_release(md_env, force_env, root_section)
150 136 : CALL update_motion_release(motion_section)
151 408 : DO ikey = 1, SIZE(keys)
152 272 : log_unit = cp_logger_get_default_io_unit(logger)
153 272 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
154 136 : motion_section, keys(ikey)), cp_p_file)) THEN
155 : ires = cp_print_key_unit_nr(logger, motion_section, TRIM(keys(ikey)), &
156 : extension=".restart.bin", &
157 : file_action="READWRITE", &
158 : file_form="UNFORMATTED", &
159 : file_position="REWIND", &
160 : file_status="UNKNOWN", &
161 272 : do_backup=(ikey == 2))
162 272 : CALL write_binary_restart(ires, log_unit, root_section, md_env, force_env)
163 : CALL cp_print_key_finished_output(ires, logger, motion_section, &
164 272 : TRIM(keys(ikey)))
165 : END IF
166 : END DO
167 : END IF
168 :
169 : CALL update_input(md_env, force_env, root_section, coords, vels, pint_env, helium_env, &
170 : save_mem=save_mem, &
171 14046 : write_binary_restart_file=write_binary_restart_file)
172 :
173 42138 : DO ikey = 1, SIZE(keys)
174 28092 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
175 14046 : motion_section, keys(ikey)), cp_p_file)) THEN
176 : ires = cp_print_key_unit_nr(logger, motion_section, TRIM(keys(ikey)), &
177 : extension=".restart", &
178 : file_position="REWIND", &
179 15526 : do_backup=(ikey == 2))
180 15526 : IF (ires > 0) THEN
181 8116 : CALL write_restart_header(ires)
182 8116 : CALL section_vals_write(root_section, unit_nr=ires, hide_root=.TRUE.)
183 : END IF
184 15526 : CALL cp_print_key_finished_output(ires, logger, motion_section, TRIM(keys(ikey)))
185 : END IF
186 : END DO
187 :
188 14046 : IF (save_mem) THEN
189 84 : CALL update_subsys_release(md_env, force_env, root_section)
190 84 : CALL update_motion_release(motion_section)
191 : END IF
192 :
193 : END IF
194 :
195 51918 : CALL timestop(handle)
196 :
197 51918 : END SUBROUTINE write_restart
198 :
199 : ! **************************************************************************************************
200 : !> \brief deallocate some sub_sections of the section subsys to save some memory
201 : !> \param md_env ...
202 : !> \param force_env ...
203 : !> \param root_section ...
204 : !> \par History
205 : !> 06.2007 created [MI]
206 : !> \author MI
207 : ! **************************************************************************************************
208 220 : SUBROUTINE update_subsys_release(md_env, force_env, root_section)
209 :
210 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
211 : TYPE(force_env_type), OPTIONAL, POINTER :: force_env
212 : TYPE(section_vals_type), POINTER :: root_section
213 :
214 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_subsys_release'
215 :
216 : CHARACTER(LEN=default_string_length) :: unit_str
217 : INTEGER :: handle, iforce_eval, myid, nforce_eval
218 220 : INTEGER, DIMENSION(:), POINTER :: i_force_eval
219 : LOGICAL :: explicit, scale, skip_vel_section
220 : TYPE(cp_subsys_type), POINTER :: subsys
221 : TYPE(force_env_type), POINTER :: my_force_b, my_force_env
222 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
223 : shell_particles
224 : TYPE(section_vals_type), POINTER :: force_env_sections, subsys_section, &
225 : work_section
226 :
227 220 : CALL timeset(routineN, handle)
228 :
229 220 : NULLIFY (core_particles, my_force_env, my_force_b, particles, &
230 220 : shell_particles, subsys, work_section)
231 :
232 220 : IF (PRESENT(md_env)) THEN
233 148 : CALL get_md_env(md_env=md_env, force_env=my_force_env)
234 72 : ELSE IF (PRESENT(force_env)) THEN
235 72 : my_force_env => force_env
236 : END IF
237 :
238 220 : IF (ASSOCIATED(my_force_env)) THEN
239 220 : NULLIFY (subsys_section)
240 220 : CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=myid)
241 : skip_vel_section = ( &
242 : (myid /= mol_dyn_run) .AND. &
243 : (myid /= mon_car_run) .AND. &
244 220 : (myid /= pint_run))
245 :
246 220 : force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
247 220 : CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval)
248 :
249 440 : DO iforce_eval = 1, nforce_eval
250 : subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", &
251 220 : i_rep_section=i_force_eval(iforce_eval))
252 220 : CALL section_vals_get(subsys_section, explicit=explicit)
253 220 : IF (.NOT. explicit) CYCLE ! Nothing to update...
254 :
255 220 : my_force_b => my_force_env
256 220 : IF (iforce_eval > 1) my_force_b => my_force_env%sub_force_env(iforce_eval - 1)%force_env
257 :
258 220 : CALL force_env_get(my_force_b, subsys=subsys)
259 :
260 : CALL cp_subsys_get(subsys, particles=particles, shell_particles=shell_particles, &
261 220 : core_particles=core_particles)
262 :
263 220 : work_section => section_vals_get_subs_vals(subsys_section, "COORD")
264 220 : CALL section_vals_get(work_section, explicit=explicit)
265 220 : IF (explicit) THEN
266 220 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
267 220 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
268 : END IF
269 220 : CALL section_vals_remove_values(work_section)
270 220 : IF (explicit) THEN
271 220 : CALL section_vals_val_set(work_section, "UNIT", c_val=unit_str)
272 220 : CALL section_vals_val_set(work_section, "SCALED", l_val=scale)
273 : END IF
274 :
275 220 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
276 220 : IF (.NOT. skip_vel_section) THEN
277 148 : CALL section_vals_remove_values(work_section)
278 : END IF
279 :
280 220 : IF (ASSOCIATED(shell_particles)) THEN
281 68 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_COORD")
282 68 : CALL section_vals_get(work_section, explicit=explicit)
283 68 : IF (explicit) THEN
284 20 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
285 20 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
286 : END IF
287 68 : CALL section_vals_remove_values(work_section)
288 68 : IF (explicit) THEN
289 20 : CALL section_vals_val_set(work_section, "UNIT", c_val=unit_str)
290 20 : CALL section_vals_val_set(work_section, "SCALED", l_val=scale)
291 : END IF
292 :
293 68 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
294 68 : IF (.NOT. skip_vel_section) THEN
295 68 : CALL section_vals_remove_values(work_section)
296 : END IF
297 : END IF
298 :
299 880 : IF (ASSOCIATED(core_particles)) THEN
300 68 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_COORD")
301 68 : CALL section_vals_get(work_section, explicit=explicit)
302 68 : IF (explicit) THEN
303 20 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
304 20 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
305 : END IF
306 68 : CALL section_vals_remove_values(work_section)
307 68 : IF (explicit) THEN
308 20 : CALL section_vals_val_set(work_section, "UNIT", c_val=unit_str)
309 20 : CALL section_vals_val_set(work_section, "SCALED", l_val=scale)
310 : END IF
311 :
312 68 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
313 68 : IF (.NOT. skip_vel_section) THEN
314 68 : CALL section_vals_remove_values(work_section)
315 : END IF
316 : END IF
317 :
318 : END DO
319 :
320 220 : DEALLOCATE (i_force_eval)
321 :
322 : END IF
323 :
324 220 : CALL timestop(handle)
325 :
326 220 : END SUBROUTINE update_subsys_release
327 :
328 : ! **************************************************************************************************
329 : !> \brief deallocate the nose subsections (coord, vel, force, mass) in the md section
330 : !> \param motion_section ...
331 : !> \par History
332 : !> 08.2007 created [MI]
333 : !> \author MI
334 : ! **************************************************************************************************
335 220 : SUBROUTINE update_motion_release(motion_section)
336 :
337 : TYPE(section_vals_type), POINTER :: motion_section
338 :
339 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_motion_release'
340 :
341 : INTEGER :: handle
342 : TYPE(section_vals_type), POINTER :: work_section
343 :
344 220 : CALL timeset(routineN, handle)
345 :
346 220 : NULLIFY (work_section)
347 :
348 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%AVERAGES%RESTART_AVERAGES")
349 220 : CALL section_vals_remove_values(work_section)
350 :
351 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%COORD")
352 220 : CALL section_vals_remove_values(work_section)
353 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%VELOCITY")
354 220 : CALL section_vals_remove_values(work_section)
355 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%MASS")
356 220 : CALL section_vals_remove_values(work_section)
357 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%FORCE")
358 220 : CALL section_vals_remove_values(work_section)
359 :
360 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%COORD")
361 220 : CALL section_vals_remove_values(work_section)
362 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%VELOCITY")
363 220 : CALL section_vals_remove_values(work_section)
364 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%MASS")
365 220 : CALL section_vals_remove_values(work_section)
366 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%FORCE")
367 220 : CALL section_vals_remove_values(work_section)
368 :
369 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%COORD")
370 220 : CALL section_vals_remove_values(work_section)
371 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%VELOCITY")
372 220 : CALL section_vals_remove_values(work_section)
373 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%MASS")
374 220 : CALL section_vals_remove_values(work_section)
375 220 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%FORCE")
376 220 : CALL section_vals_remove_values(work_section)
377 :
378 220 : CALL timestop(handle)
379 :
380 220 : END SUBROUTINE update_motion_release
381 :
382 : ! **************************************************************************************************
383 : !> \brief Updates the whole input file for the restart
384 : !> \param md_env ...
385 : !> \param force_env ...
386 : !> \param root_section ...
387 : !> \param coords ...
388 : !> \param vels ...
389 : !> \param pint_env ...
390 : !> \param helium_env ...
391 : !> \param save_mem ...
392 : !> \param write_binary_restart_file ...
393 : !> \par History
394 : !> 01.2006 created [teo]
395 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
396 : !> \author Teodoro Laino
397 : ! **************************************************************************************************
398 14046 : SUBROUTINE update_input(md_env, force_env, root_section, coords, vels, pint_env, &
399 : helium_env, save_mem, write_binary_restart_file)
400 :
401 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
402 : TYPE(force_env_type), OPTIONAL, POINTER :: force_env
403 : TYPE(section_vals_type), POINTER :: root_section
404 : TYPE(neb_var_type), OPTIONAL, POINTER :: coords, vels
405 : TYPE(pint_env_type), INTENT(IN), OPTIONAL :: pint_env
406 : TYPE(helium_solvent_p_type), DIMENSION(:), &
407 : OPTIONAL, POINTER :: helium_env
408 : LOGICAL, INTENT(IN), OPTIONAL :: save_mem, write_binary_restart_file
409 :
410 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_input'
411 :
412 : INTEGER :: handle
413 : LOGICAL :: do_respa, lcond, my_save_mem, &
414 : my_write_binary_restart_file
415 : TYPE(cp_logger_type), POINTER :: logger
416 : TYPE(force_env_type), POINTER :: my_force_env
417 : TYPE(section_vals_type), POINTER :: motion_section
418 : TYPE(simpar_type), POINTER :: simpar
419 :
420 14046 : CALL timeset(routineN, handle)
421 :
422 14046 : NULLIFY (logger, motion_section, my_force_env)
423 :
424 : IF (PRESENT(save_mem)) THEN
425 : my_save_mem = save_mem
426 : ELSE
427 : my_save_mem = .FALSE.
428 : END IF
429 :
430 14046 : IF (PRESENT(write_binary_restart_file)) THEN
431 14046 : my_write_binary_restart_file = write_binary_restart_file
432 : ELSE
433 0 : my_write_binary_restart_file = .FALSE.
434 : END IF
435 :
436 14046 : logger => cp_get_default_logger()
437 :
438 : ! Can handle md_env or force_env
439 14046 : lcond = PRESENT(md_env) .OR. PRESENT(force_env) .OR. PRESENT(pint_env) .OR. PRESENT(helium_env)
440 : IF (lcond) THEN
441 13916 : IF (PRESENT(md_env)) THEN
442 5526 : CALL get_md_env(md_env=md_env, force_env=my_force_env)
443 8390 : ELSE IF (PRESENT(force_env)) THEN
444 7752 : my_force_env => force_env
445 : END IF
446 : ! The real restart setting...
447 13916 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
448 : CALL update_motion(motion_section, &
449 : md_env=md_env, &
450 : force_env=my_force_env, &
451 : logger=logger, &
452 : coords=coords, &
453 : vels=vels, &
454 : pint_env=pint_env, &
455 : helium_env=helium_env, &
456 : save_mem=my_save_mem, &
457 13916 : write_binary_restart_file=my_write_binary_restart_file)
458 : ! Update one force_env_section per time..
459 13916 : IF (ASSOCIATED(my_force_env)) THEN
460 13278 : do_respa = .FALSE.
461 : ! Do respa only in case of RESPA MD
462 13278 : IF (PRESENT(md_env)) THEN
463 5526 : CALL get_md_env(md_env=md_env, simpar=simpar)
464 5526 : IF (simpar%do_respa) THEN
465 6 : do_respa = .TRUE.
466 : END IF
467 : END IF
468 :
469 : CALL update_force_eval(force_env=my_force_env, &
470 : root_section=root_section, &
471 : write_binary_restart_file=my_write_binary_restart_file, &
472 13278 : respa=do_respa)
473 :
474 : END IF
475 : END IF
476 :
477 14046 : CALL timestop(handle)
478 :
479 14046 : END SUBROUTINE update_input
480 :
481 : ! **************************************************************************************************
482 : !> \brief Updates the motion section of the input file
483 : !> \param motion_section ...
484 : !> \param md_env ...
485 : !> \param force_env ...
486 : !> \param logger ...
487 : !> \param coords ...
488 : !> \param vels ...
489 : !> \param pint_env ...
490 : !> \param helium_env ...
491 : !> \param save_mem ...
492 : !> \param write_binary_restart_file ...
493 : !> \par History
494 : !> 01.2006 created [teo]
495 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
496 : !> \author Teodoro Laino
497 : ! **************************************************************************************************
498 125244 : SUBROUTINE update_motion(motion_section, md_env, force_env, logger, &
499 : coords, vels, pint_env, helium_env, save_mem, &
500 : write_binary_restart_file)
501 :
502 : TYPE(section_vals_type), POINTER :: motion_section
503 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
504 : TYPE(force_env_type), POINTER :: force_env
505 : TYPE(cp_logger_type), POINTER :: logger
506 : TYPE(neb_var_type), OPTIONAL, POINTER :: coords, vels
507 : TYPE(pint_env_type), INTENT(IN), OPTIONAL :: pint_env
508 : TYPE(helium_solvent_p_type), DIMENSION(:), &
509 : OPTIONAL, POINTER :: helium_env
510 : LOGICAL, INTENT(IN), OPTIONAL :: save_mem, write_binary_restart_file
511 :
512 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_motion'
513 :
514 : INTEGER :: counter, handle, handle2, i, irep, isec, &
515 : j, nhc_len, tot_nhcneed
516 13916 : INTEGER, DIMENSION(:), POINTER :: walkers_status
517 : INTEGER, POINTER :: itimes
518 : LOGICAL :: my_save_mem, my_write_binary_restart_file
519 13916 : REAL(KIND=dp), DIMENSION(:), POINTER :: buffer, eta, fnhc, mnhc, veta, wrk
520 : REAL(KIND=dp), POINTER :: constant, t
521 : TYPE(average_quantities_type), POINTER :: averages
522 : TYPE(cp_subsys_type), POINTER :: subsys
523 : TYPE(lnhc_parameters_type), POINTER :: nhc
524 : TYPE(meta_env_type), POINTER :: meta_env
525 : TYPE(mp_para_env_type), POINTER :: para_env
526 13916 : TYPE(npt_info_type), POINTER :: npt(:, :)
527 : TYPE(particle_list_type), POINTER :: particles
528 : TYPE(section_vals_type), POINTER :: replica_section, work_section
529 : TYPE(simpar_type), POINTER :: simpar
530 : TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_part, &
531 : thermostat_shell
532 :
533 13916 : CALL timeset(routineN, handle)
534 13916 : NULLIFY (logger, thermostat_part, thermostat_baro, npt, para_env, nhc, &
535 13916 : work_section, thermostat_shell, t, averages, constant, &
536 13916 : walkers_status, itimes, meta_env, simpar)
537 13916 : NULLIFY (particles)
538 13916 : NULLIFY (subsys)
539 13916 : IF (PRESENT(md_env)) THEN
540 : CALL get_md_env(md_env=md_env, &
541 : thermostat_part=thermostat_part, &
542 : thermostat_baro=thermostat_baro, &
543 : thermostat_shell=thermostat_shell, &
544 : npt=npt, &
545 : t=t, &
546 : constant=constant, &
547 : itimes=itimes, &
548 : simpar=simpar, &
549 : averages=averages, &
550 5526 : para_env=para_env)
551 : ELSE
552 8390 : IF (ASSOCIATED(force_env)) THEN
553 7752 : para_env => force_env%para_env
554 638 : ELSE IF (PRESENT(pint_env)) THEN
555 600 : para_env => pint_env%logger%para_env
556 38 : ELSE IF (PRESENT(helium_env)) THEN
557 : ! Only needed in case that pure helium is simulated
558 : ! In this case write_restart is called only by processors
559 : ! with associated helium_env
560 38 : para_env => helium_env(1)%helium%logger%para_env
561 : ELSE
562 0 : CPABORT("No valid para_env present")
563 : END IF
564 : END IF
565 :
566 13916 : IF (ASSOCIATED(force_env)) THEN
567 13278 : meta_env => force_env%meta_env
568 : END IF
569 :
570 : IF (PRESENT(save_mem)) THEN
571 13916 : my_save_mem = save_mem
572 : ELSE
573 13916 : my_save_mem = .FALSE.
574 : END IF
575 :
576 13916 : IF (PRESENT(write_binary_restart_file)) THEN
577 13916 : my_write_binary_restart_file = write_binary_restart_file
578 : ELSE
579 : my_write_binary_restart_file = .FALSE.
580 : END IF
581 :
582 13916 : CALL timeset(routineN//"_COUNTERS", handle2)
583 13916 : IF (ASSOCIATED(itimes)) THEN
584 5526 : IF (itimes >= 0) THEN
585 5526 : CALL section_vals_val_set(motion_section, "MD%STEP_START_VAL", i_val=itimes)
586 5526 : CPASSERT(ASSOCIATED(t))
587 5526 : CALL section_vals_val_set(motion_section, "MD%TIME_START_VAL", r_val=t)
588 : END IF
589 : END IF
590 13916 : IF (ASSOCIATED(constant)) THEN
591 5526 : CALL section_vals_val_set(motion_section, "MD%ECONS_START_VAL", r_val=constant)
592 : END IF
593 13916 : CALL timestop(handle2)
594 : ! AVERAGES
595 13916 : CALL timeset(routineN//"_AVERAGES", handle2)
596 13916 : IF (ASSOCIATED(averages)) THEN
597 5526 : IF ((averages%do_averages) .AND. (averages%itimes_start /= -1)) THEN
598 5518 : work_section => section_vals_get_subs_vals(motion_section, "MD%AVERAGES")
599 5518 : CALL section_vals_val_set(work_section, "_SECTION_PARAMETERS_", l_val=averages%do_averages)
600 5518 : work_section => section_vals_get_subs_vals(motion_section, "MD%AVERAGES%RESTART_AVERAGES")
601 5518 : CALL section_vals_val_set(work_section, "ITIMES_START", i_val=averages%itimes_start)
602 5518 : CALL section_vals_val_set(work_section, "AVECPU", r_val=averages%avecpu)
603 5518 : CALL section_vals_val_set(work_section, "AVEHUGONIOT", r_val=averages%avehugoniot)
604 5518 : CALL section_vals_val_set(work_section, "AVETEMP_BARO", r_val=averages%avetemp_baro)
605 5518 : CALL section_vals_val_set(work_section, "AVEPOT", r_val=averages%avepot)
606 5518 : CALL section_vals_val_set(work_section, "AVEKIN", r_val=averages%avekin)
607 5518 : CALL section_vals_val_set(work_section, "AVETEMP", r_val=averages%avetemp)
608 5518 : CALL section_vals_val_set(work_section, "AVEKIN_QM", r_val=averages%avekin_qm)
609 5518 : CALL section_vals_val_set(work_section, "AVETEMP_QM", r_val=averages%avetemp_qm)
610 5518 : CALL section_vals_val_set(work_section, "AVEVOL", r_val=averages%avevol)
611 5518 : CALL section_vals_val_set(work_section, "AVECELL_A", r_val=averages%aveca)
612 5518 : CALL section_vals_val_set(work_section, "AVECELL_B", r_val=averages%avecb)
613 5518 : CALL section_vals_val_set(work_section, "AVECELL_C", r_val=averages%avecc)
614 5518 : CALL section_vals_val_set(work_section, "AVEALPHA", r_val=averages%aveal)
615 5518 : CALL section_vals_val_set(work_section, "AVEBETA", r_val=averages%avebe)
616 5518 : CALL section_vals_val_set(work_section, "AVEGAMMA", r_val=averages%avega)
617 5518 : CALL section_vals_val_set(work_section, "AVE_ECONS", r_val=averages%econs)
618 5518 : CALL section_vals_val_set(work_section, "AVE_PRESS", r_val=averages%avepress)
619 5518 : CALL section_vals_val_set(work_section, "AVE_PXX", r_val=averages%avepxx)
620 : ! Virial averages
621 5518 : IF (ASSOCIATED(averages%virial)) THEN
622 0 : ALLOCATE (buffer(9))
623 0 : buffer = RESHAPE(averages%virial%pv_total, [9])
624 0 : CALL section_vals_val_set(work_section, "AVE_PV_TOT", r_vals_ptr=buffer)
625 :
626 0 : ALLOCATE (buffer(9))
627 0 : buffer = RESHAPE(averages%virial%pv_virial, [9])
628 0 : CALL section_vals_val_set(work_section, "AVE_PV_VIR", r_vals_ptr=buffer)
629 :
630 0 : ALLOCATE (buffer(9))
631 0 : buffer = RESHAPE(averages%virial%pv_kinetic, [9])
632 0 : CALL section_vals_val_set(work_section, "AVE_PV_KIN", r_vals_ptr=buffer)
633 :
634 0 : ALLOCATE (buffer(9))
635 0 : buffer = RESHAPE(averages%virial%pv_constraint, [9])
636 0 : CALL section_vals_val_set(work_section, "AVE_PV_CNSTR", r_vals_ptr=buffer)
637 :
638 0 : ALLOCATE (buffer(9))
639 0 : buffer = RESHAPE(averages%virial%pv_xc, [9])
640 0 : CALL section_vals_val_set(work_section, "AVE_PV_XC", r_vals_ptr=buffer)
641 :
642 0 : ALLOCATE (buffer(9))
643 0 : buffer = RESHAPE(averages%virial%pv_fock_4c, [9])
644 0 : CALL section_vals_val_set(work_section, "AVE_PV_FOCK_4C", r_vals_ptr=buffer)
645 : END IF
646 : ! Colvars averages
647 5518 : IF (SIZE(averages%avecolvar) > 0) THEN
648 6 : ALLOCATE (buffer(SIZE(averages%avecolvar)))
649 196 : buffer = averages%avecolvar
650 2 : CALL section_vals_val_set(work_section, "AVE_COLVARS", r_vals_ptr=buffer)
651 : END IF
652 5518 : IF (SIZE(averages%aveMmatrix) > 0) THEN
653 6 : ALLOCATE (buffer(SIZE(averages%aveMmatrix)))
654 9220 : buffer = averages%aveMmatrix
655 2 : CALL section_vals_val_set(work_section, "AVE_MMATRIX", r_vals_ptr=buffer)
656 : END IF
657 : END IF
658 : END IF
659 13916 : CALL timestop(handle2)
660 :
661 : ! SAVE THERMOSTAT target TEMPERATURE when doing TEMPERATURE_ANNEALING
662 13916 : IF (PRESENT(md_env)) THEN
663 5526 : IF (ASSOCIATED(simpar)) THEN
664 5526 : IF (simpar%temperature_annealing .AND. ABS(1._dp - simpar%f_temperature_annealing) > 1.E-10_dp) THEN
665 4 : CALL section_vals_val_set(motion_section, "MD%TEMPERATURE", r_val=simpar%temp_ext)
666 : END IF
667 : END IF
668 : END IF
669 :
670 : ! PARTICLE THERMOSTAT
671 13916 : CALL timeset(routineN//"_THERMOSTAT_PARTICLES", handle2)
672 13916 : IF (ASSOCIATED(thermostat_part)) THEN
673 1042 : IF (thermostat_part%type_of_thermostat == do_thermo_nose) THEN
674 : ! Restart of Nose-Hoover Thermostat for Particles
675 728 : IF (.NOT. my_write_binary_restart_file) THEN
676 640 : nhc => thermostat_part%nhc
677 640 : CALL collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
678 640 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE")
679 640 : CALL set_template_restart(work_section, eta, veta, fnhc, mnhc)
680 : END IF
681 314 : ELSE IF (thermostat_part%type_of_thermostat == do_thermo_csvr) THEN
682 : ! Restart of CSVR Thermostat for Particles
683 294 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%CSVR")
684 294 : CALL dump_csvr_restart_info(thermostat_part%csvr, para_env, work_section)
685 20 : ELSE IF (thermostat_part%type_of_thermostat == do_thermo_al) THEN
686 : ! Restart of AD_LANGEVIN Thermostat for Particles
687 0 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%AD_LANGEVIN")
688 0 : CALL dump_al_restart_info(thermostat_part%al, para_env, work_section)
689 20 : ELSE IF (thermostat_part%type_of_thermostat == do_thermo_gle) THEN
690 : ! Restart of GLE Thermostat for Particles
691 20 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%GLE")
692 20 : CALL dump_gle_restart_info(thermostat_part%gle, para_env, work_section)
693 : END IF
694 : END IF
695 13916 : CALL timestop(handle2)
696 :
697 : ! BAROSTAT - THERMOSTAT
698 13916 : CALL timeset(routineN//"_BAROSTAT", handle2)
699 13916 : IF (ASSOCIATED(thermostat_baro)) THEN
700 336 : IF (thermostat_baro%type_of_thermostat == do_thermo_nose) THEN
701 : ! Restart of Nose-Hoover Thermostat for Barostat
702 252 : nhc => thermostat_baro%nhc
703 252 : nhc_len = SIZE(nhc%nvt, 1)
704 252 : tot_nhcneed = nhc%glob_num_nhc
705 756 : ALLOCATE (eta(tot_nhcneed*nhc_len))
706 504 : ALLOCATE (veta(tot_nhcneed*nhc_len))
707 504 : ALLOCATE (fnhc(tot_nhcneed*nhc_len))
708 504 : ALLOCATE (mnhc(tot_nhcneed*nhc_len))
709 252 : counter = 0
710 1148 : DO i = 1, SIZE(nhc%nvt, 1)
711 2044 : DO j = 1, SIZE(nhc%nvt, 2)
712 896 : counter = counter + 1
713 896 : eta(counter) = nhc%nvt(i, j)%eta
714 896 : veta(counter) = nhc%nvt(i, j)%v
715 896 : fnhc(counter) = nhc%nvt(i, j)%f
716 1792 : mnhc(counter) = nhc%nvt(i, j)%mass
717 : END DO
718 : END DO
719 252 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE")
720 252 : CALL set_template_restart(work_section, eta, veta, fnhc, mnhc)
721 84 : ELSE IF (thermostat_baro%type_of_thermostat == do_thermo_csvr) THEN
722 : ! Restart of CSVR Thermostat for Barostat
723 84 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%CSVR")
724 84 : CALL dump_csvr_restart_info(thermostat_baro%csvr, para_env, work_section)
725 : END IF
726 : END IF
727 13916 : CALL timestop(handle2)
728 :
729 : ! BAROSTAT
730 13916 : CALL timeset(routineN//"_NPT", handle2)
731 13916 : IF (ASSOCIATED(npt)) THEN
732 1194 : ALLOCATE (veta(SIZE(npt, 1)*SIZE(npt, 2)))
733 796 : ALLOCATE (mnhc(SIZE(npt, 1)*SIZE(npt, 2)))
734 398 : counter = 0
735 1124 : DO i = 1, SIZE(npt, 1)
736 2834 : DO j = 1, SIZE(npt, 2)
737 1710 : counter = counter + 1
738 1710 : veta(counter) = npt(i, j)%v
739 2436 : mnhc(counter) = npt(i, j)%mass
740 : END DO
741 : END DO
742 398 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT")
743 398 : CALL set_template_restart(work_section, veta=veta, mnhc=mnhc)
744 : END IF
745 13916 : CALL timestop(handle2)
746 :
747 : ! SHELL THERMOSTAT
748 13916 : CALL timeset(routineN//"_THERMOSTAT_SHELL", handle2)
749 13916 : IF (ASSOCIATED(thermostat_shell)) THEN
750 160 : IF (thermostat_shell%type_of_thermostat == do_thermo_nose) THEN
751 : ! Restart of Nose-Hoover Thermostat for Shell Particles
752 136 : IF (.NOT. my_write_binary_restart_file) THEN
753 124 : nhc => thermostat_shell%nhc
754 124 : CALL collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
755 124 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE")
756 124 : CALL set_template_restart(work_section, eta, veta, fnhc, mnhc)
757 : END IF
758 24 : ELSE IF (thermostat_shell%type_of_thermostat == do_thermo_csvr) THEN
759 24 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%CSVR")
760 : ! Restart of CSVR Thermostat for Shell Particles
761 24 : CALL dump_csvr_restart_info(thermostat_shell%csvr, para_env, work_section)
762 : END IF
763 : END IF
764 13916 : CALL timestop(handle2)
765 :
766 13916 : CALL timeset(routineN//"_META", handle2)
767 13916 : IF (ASSOCIATED(meta_env)) THEN
768 : CALL section_vals_val_set(meta_env%metadyn_section, "STEP_START_VAL", &
769 990 : i_val=meta_env%n_steps)
770 : CALL section_vals_val_set(meta_env%metadyn_section, "NHILLS_START_VAL", &
771 990 : i_val=meta_env%hills_env%n_hills)
772 : !RG Adaptive hills
773 : CALL section_vals_val_set(meta_env%metadyn_section, "MIN_DISP", &
774 990 : r_val=meta_env%hills_env%min_disp)
775 : CALL section_vals_val_set(meta_env%metadyn_section, "OLD_HILL_NUMBER", &
776 990 : i_val=meta_env%hills_env%old_hill_number)
777 : CALL section_vals_val_set(meta_env%metadyn_section, "OLD_HILL_STEP", &
778 990 : i_val=meta_env%hills_env%old_hill_step)
779 : !RG Adaptive hills
780 990 : IF (meta_env%do_hills .AND. meta_env%hills_env%n_hills /= 0) THEN
781 784 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_POS")
782 784 : CALL meta_hills_val_set_ss(work_section, meta_env)
783 784 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_SCALE")
784 784 : CALL meta_hills_val_set_ds(work_section, meta_env)
785 784 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_HEIGHT")
786 784 : CALL meta_hills_val_set_ww(work_section, meta_env)
787 784 : IF (meta_env%well_tempered) THEN
788 2 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_INVDT")
789 2 : CALL meta_hills_val_set_dt(work_section, meta_env)
790 : END IF
791 : END IF
792 990 : IF (meta_env%extended_lagrange) THEN
793 : CALL section_vals_val_set(meta_env%metadyn_section, "COLVAR_AVG_TEMPERATURE_RESTART", &
794 130 : r_val=meta_env%avg_temp)
795 130 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
796 290 : DO irep = 1, meta_env%n_colvar
797 : CALL section_vals_val_set(work_section, "_DEFAULT_KEYWORD_", r_val=meta_env%metavar(irep)%ss0, &
798 290 : i_rep_val=irep)
799 : END DO
800 130 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
801 290 : DO irep = 1, meta_env%n_colvar
802 : CALL section_vals_val_set(work_section, "_DEFAULT_KEYWORD_", r_val=meta_env%metavar(irep)%vvp, &
803 290 : i_rep_val=irep)
804 : END DO
805 :
806 130 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS")
807 290 : DO irep = 1, meta_env%n_colvar
808 : CALL section_vals_val_set(work_section, "_DEFAULT_KEYWORD_", r_val=meta_env%metavar(irep)%ss, &
809 290 : i_rep_val=irep)
810 : END DO
811 130 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_FS")
812 290 : DO irep = 1, meta_env%n_colvar
813 : CALL section_vals_val_set(work_section, "_DEFAULT_KEYWORD_", r_val=meta_env%metavar(irep)%ff_s, &
814 290 : i_rep_val=irep)
815 : END DO
816 :
817 : END IF
818 : ! Multiple Walkers
819 990 : IF (meta_env%do_multiple_walkers) THEN
820 636 : ALLOCATE (walkers_status(meta_env%multiple_walkers%walkers_tot_nr))
821 1272 : walkers_status = meta_env%multiple_walkers%walkers_status
822 212 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "MULTIPLE_WALKERS")
823 212 : CALL section_vals_val_set(work_section, "WALKERS_STATUS", i_vals_ptr=walkers_status)
824 : END IF
825 : END IF
826 13916 : CALL timestop(handle2)
827 13916 : CALL timeset(routineN//"_NEB", handle2)
828 13916 : IF (PRESENT(coords) .OR. (PRESENT(vels))) THEN
829 : ! Update NEB section
830 578 : replica_section => section_vals_get_subs_vals(motion_section, "BAND%REPLICA")
831 578 : CALL force_env_get(force_env, subsys=subsys)
832 578 : CALL cp_subsys_get(subsys, particles=particles)
833 578 : IF (PRESENT(coords)) THEN
834 : ! Allocate possible missing sections
835 52 : DO
836 630 : IF (coords%size_wrk(2) <= SIZE(replica_section%values, 2)) EXIT
837 52 : CALL section_vals_add_values(replica_section)
838 : END DO
839 : ! Write Values
840 4152 : DO isec = 1, coords%size_wrk(2)
841 3574 : CALL section_vals_val_unset(replica_section, "COORD_FILE_NAME", i_rep_section=isec)
842 3574 : work_section => section_vals_get_subs_vals3(replica_section, "COORD", i_rep_section=isec)
843 : CALL section_neb_coord_val_set(work_section, coords%xyz(:, isec), SIZE(coords%xyz, 1), 3*SIZE(particles%els), &
844 3574 : 3, particles%els, angstrom)
845 : ! Update Collective Variables
846 4152 : IF (coords%in_use == do_band_collective) THEN
847 360 : ALLOCATE (wrk(coords%size_wrk(1)))
848 480 : wrk = coords%wrk(:, isec)
849 : CALL section_vals_val_set(replica_section, "COLLECTIVE", r_vals_ptr=wrk, &
850 120 : i_rep_section=isec)
851 : END IF
852 : END DO
853 : END IF
854 578 : IF (PRESENT(vels)) THEN
855 578 : CALL force_env_get(force_env, subsys=subsys)
856 578 : CALL cp_subsys_get(subsys, particles=particles)
857 : ! Allocate possible missing sections
858 0 : DO
859 578 : IF (vels%size_wrk(2) <= SIZE(replica_section%values, 2)) EXIT
860 0 : CALL section_vals_add_values(replica_section)
861 : END DO
862 : ! Write Values
863 4152 : DO isec = 1, vels%size_wrk(2)
864 3574 : work_section => section_vals_get_subs_vals3(replica_section, "VELOCITY", i_rep_section=isec)
865 4152 : IF (vels%in_use == do_band_collective) THEN
866 : CALL section_neb_coord_val_set(work_section, vels%wrk(:, isec), SIZE(vels%wrk, 1), SIZE(vels%wrk, 1), &
867 120 : 1, particles%els, 1.0_dp)
868 : ELSE
869 : CALL section_neb_coord_val_set(work_section, vels%wrk(:, isec), SIZE(vels%wrk, 1), 3*SIZE(particles%els), &
870 3454 : 3, particles%els, 1.0_dp)
871 : END IF
872 : END DO
873 : END IF
874 : END IF
875 13916 : CALL timestop(handle2)
876 :
877 13916 : IF (PRESENT(pint_env)) THEN
878 : ! Update PINT section
879 600 : CALL update_motion_pint(motion_section, pint_env)
880 : END IF
881 :
882 13916 : IF (PRESENT(helium_env)) THEN
883 : ! Update HELIUM section
884 110 : CALL update_motion_helium(helium_env)
885 : END IF
886 :
887 13916 : CALL timestop(handle)
888 :
889 13916 : END SUBROUTINE update_motion
890 :
891 : ! ***************************************************************************
892 : !> \brief Update PINT section in the input structure
893 : !> \param motion_section ...
894 : !> \param pint_env ...
895 : !> \date 2010-10-13
896 : !> \author Lukasz Walewski <Lukasz.Walewski@ruhr-uni-bochum.de>
897 : ! **************************************************************************************************
898 600 : SUBROUTINE update_motion_pint(motion_section, pint_env)
899 :
900 : TYPE(section_vals_type), POINTER :: motion_section
901 : TYPE(pint_env_type), INTENT(IN) :: pint_env
902 :
903 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_motion_pint'
904 :
905 : CHARACTER(LEN=rng_record_length) :: rng_record
906 : INTEGER :: handle, i, iatom, ibead, inos, isp
907 : INTEGER, DIMENSION(rng_record_length, 1) :: ascii
908 : LOGICAL :: explicit
909 600 : REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals
910 : TYPE(section_vals_type), POINTER :: pint_section, tmpsec
911 :
912 600 : CALL timeset(routineN, handle)
913 :
914 600 : pint_section => section_vals_get_subs_vals(motion_section, "PINT")
915 600 : CALL section_vals_val_set(pint_section, "ITERATION", i_val=pint_env%iter)
916 :
917 : ! allocate memory for COORDs and VELOCITYs if the BEADS section was not
918 : ! explicitly given in the input (this is actually done only once since
919 : ! after section_vals_add_values section becomes explicit)
920 600 : NULLIFY (tmpsec)
921 600 : tmpsec => section_vals_get_subs_vals(pint_section, "BEADS")
922 600 : CALL section_vals_get(tmpsec, explicit=explicit)
923 600 : IF (.NOT. explicit) THEN
924 46 : CALL section_vals_add_values(tmpsec)
925 : END IF
926 :
927 : ! update bead coordinates in the global input structure
928 600 : NULLIFY (r_vals)
929 1800 : ALLOCATE (r_vals(pint_env%p*pint_env%ndim))
930 :
931 600 : i = 1
932 600 : CALL pint_u2x(pint_env)
933 98160 : DO iatom = 1, pint_env%ndim
934 501072 : DO ibead = 1, pint_env%p
935 402912 : r_vals(i) = pint_env%x(ibead, iatom)
936 500472 : i = i + 1
937 : END DO
938 : END DO
939 : CALL section_vals_val_set(pint_section, "BEADS%COORD%_DEFAULT_KEYWORD_", &
940 600 : r_vals_ptr=r_vals)
941 :
942 : ! update bead velocities in the global input structure
943 600 : NULLIFY (r_vals)
944 1800 : ALLOCATE (r_vals(pint_env%p*pint_env%ndim))
945 600 : i = 1
946 600 : CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
947 98160 : DO iatom = 1, pint_env%ndim
948 501072 : DO ibead = 1, pint_env%p
949 402912 : r_vals(i) = pint_env%v(ibead, iatom)
950 500472 : i = i + 1
951 : END DO
952 : END DO
953 : CALL section_vals_val_set(pint_section, "BEADS%VELOCITY%_DEFAULT_KEYWORD_", &
954 600 : r_vals_ptr=r_vals)
955 :
956 600 : IF (pint_env%pimd_thermostat == thermostat_nose) THEN
957 :
958 : ! allocate memory for COORDs and VELOCITYs if the NOSE section was not
959 : ! explicitly given in the input (this is actually done only once since
960 : ! after section_vals_add_values section becomes explicit)
961 226 : NULLIFY (tmpsec)
962 226 : tmpsec => section_vals_get_subs_vals(pint_section, "NOSE")
963 226 : CALL section_vals_get(tmpsec, explicit=explicit)
964 226 : IF (.NOT. explicit) THEN
965 0 : CALL section_vals_add_values(tmpsec)
966 : END IF
967 :
968 : ! update thermostat coordinates in the global input structure
969 226 : NULLIFY (r_vals)
970 678 : ALLOCATE (r_vals(pint_env%p*pint_env%ndim*pint_env%nnos))
971 226 : i = 1
972 2440 : DO iatom = 1, pint_env%ndim
973 17056 : DO ibead = 1, pint_env%p
974 63558 : DO inos = 1, pint_env%nnos
975 46728 : r_vals(i) = pint_env%tx(inos, ibead, iatom)
976 61344 : i = i + 1
977 : END DO
978 : END DO
979 : END DO
980 : CALL section_vals_val_set(pint_section, "NOSE%COORD%_DEFAULT_KEYWORD_", &
981 226 : r_vals_ptr=r_vals)
982 :
983 : ! update thermostat velocities in the global input structure
984 226 : NULLIFY (r_vals)
985 678 : ALLOCATE (r_vals(pint_env%p*pint_env%ndim*pint_env%nnos))
986 226 : i = 1
987 2440 : DO iatom = 1, pint_env%ndim
988 17056 : DO ibead = 1, pint_env%p
989 63558 : DO inos = 1, pint_env%nnos
990 46728 : r_vals(i) = pint_env%tv(inos, ibead, iatom)
991 61344 : i = i + 1
992 : END DO
993 : END DO
994 : END DO
995 : CALL section_vals_val_set(pint_section, "NOSE%VELOCITY%_DEFAULT_KEYWORD_", &
996 452 : r_vals_ptr=r_vals)
997 :
998 : ELSE IF (pint_env%pimd_thermostat == thermostat_gle) THEN
999 :
1000 0 : NULLIFY (tmpsec)
1001 0 : tmpsec => section_vals_get_subs_vals(pint_section, "GLE")
1002 0 : CALL dump_gle_restart_info(pint_env%gle, pint_env%replicas%para_env, tmpsec)
1003 :
1004 : ELSE IF (pint_env%pimd_thermostat == thermostat_pile) THEN
1005 : tmpsec => section_vals_get_subs_vals(pint_section, &
1006 302 : "PILE%RNG_INIT")
1007 302 : CALL pint_env%pile_therm%gaussian_rng_stream%dump(rng_record)
1008 302 : CALL string_to_ascii(rng_record, ascii(:, 1))
1009 : CALL section_rng_val_set(rng_section=tmpsec, nsize=1, &
1010 302 : ascii=ascii)
1011 302 : tmpsec => section_vals_get_subs_vals(pint_section, "PILE")
1012 : CALL section_vals_val_set(tmpsec, "THERMOSTAT_ENERGY", &
1013 302 : r_val=pint_env%e_pile)
1014 : ELSE IF (pint_env%pimd_thermostat == thermostat_qtb) THEN
1015 : tmpsec => section_vals_get_subs_vals(pint_section, &
1016 20 : "QTB%RNG_INIT")
1017 : CALL string_to_ascii(pint_env%qtb_therm%rng_status(1), &
1018 20 : ascii(:, 1))
1019 : CALL section_rng_val_set(rng_section=tmpsec, nsize=1, &
1020 20 : ascii=ascii)
1021 20 : tmpsec => section_vals_get_subs_vals(pint_section, "QTB")
1022 : CALL section_vals_val_set(tmpsec, "THERMOSTAT_ENERGY", &
1023 20 : r_val=pint_env%e_qtb)
1024 : ELSE IF (pint_env%pimd_thermostat == thermostat_piglet) THEN
1025 : tmpsec => section_vals_get_subs_vals(pint_section, &
1026 0 : "PIGLET%RNG_INIT")
1027 0 : CALL pint_env%piglet_therm%gaussian_rng_stream%dump(rng_record)
1028 0 : CALL string_to_ascii(rng_record, ascii(:, 1))
1029 : CALL section_rng_val_set(rng_section=tmpsec, nsize=1, &
1030 0 : ascii=ascii)
1031 0 : tmpsec => section_vals_get_subs_vals(pint_section, "PIGLET")
1032 : CALL section_vals_val_set(tmpsec, "THERMOSTAT_ENERGY", &
1033 0 : r_val=pint_env%e_piglet)
1034 : ! update thermostat velocities in the global input structure
1035 0 : NULLIFY (r_vals)
1036 : ALLOCATE (r_vals((pint_env%piglet_therm%nsp1 - 1)* &
1037 : pint_env%piglet_therm%ndim* &
1038 0 : pint_env%piglet_therm%p))
1039 0 : i = 1
1040 0 : DO isp = 2, pint_env%piglet_therm%nsp1
1041 0 : DO ibead = 1, pint_env%piglet_therm%p*pint_env%piglet_therm%ndim
1042 0 : r_vals(i) = pint_env%piglet_therm%smalls(isp, ibead)
1043 0 : i = i + 1
1044 : END DO
1045 : END DO
1046 : CALL section_vals_val_set(pint_section, "PIGLET%EXTRA_DOF%_DEFAULT_KEYWORD_", &
1047 0 : r_vals_ptr=r_vals)
1048 : END IF
1049 :
1050 600 : CALL timestop(handle)
1051 :
1052 1200 : END SUBROUTINE update_motion_pint
1053 :
1054 : ! ***************************************************************************
1055 : !> \brief Update HELIUM section in the input structure.
1056 : !> \param helium_env ...
1057 : !> \date 2009-11-12
1058 : !> \parm History
1059 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
1060 : !> \author Lukasz Walewski <Lukasz.Walewski@ruhr-uni-bochum.de>
1061 : !> \note Transfer the current helium state from the runtime environment
1062 : !> to the input structure, so that it can be used for I/O, etc.
1063 : !> \note Moved from the helium_io module directly, might be done better way
1064 : ! **************************************************************************************************
1065 110 : SUBROUTINE update_motion_helium(helium_env)
1066 :
1067 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1068 :
1069 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_motion_helium'
1070 :
1071 : CHARACTER(LEN=default_string_length) :: err_str, stmp
1072 : INTEGER :: handle, i, itmp, iweight, msglen, &
1073 : nsteps, off, offset, reqlen
1074 110 : INTEGER, DIMENSION(:), POINTER :: int_msg_gather
1075 : LOGICAL :: lbf
1076 : REAL(kind=dp) :: bf, bu, invproc
1077 : REAL(kind=dp), DIMENSION(3, 2) :: bg, cg, ig
1078 110 : REAL(kind=dp), DIMENSION(:), POINTER :: real_msg, real_msg_gather
1079 : TYPE(cp_logger_type), POINTER :: logger
1080 :
1081 110 : CALL timeset(routineN, handle)
1082 :
1083 : !CPASSERT(ASSOCIATED(helium_env))
1084 :
1085 110 : NULLIFY (logger)
1086 110 : logger => cp_get_default_logger()
1087 :
1088 110 : IF (ASSOCIATED(helium_env)) THEN
1089 : ! determine offset for arrays
1090 105 : offset = 0
1091 155 : DO i = 1, logger%para_env%mepos
1092 155 : offset = offset + helium_env(1)%env_all(i)
1093 : END DO
1094 :
1095 105 : IF (.NOT. helium_env(1)%helium%solute_present) THEN
1096 : ! update iteration number
1097 38 : itmp = logger%iter_info%iteration(2)
1098 : CALL section_vals_val_set( &
1099 : helium_env(1)%helium%input, &
1100 : "MOTION%PINT%ITERATION", &
1101 38 : i_val=itmp)
1102 : ! else - PINT will do that
1103 : END IF
1104 :
1105 : !
1106 : ! save coordinates
1107 : !
1108 : ! allocate the buffer to be passed and fill it with local coords at each
1109 : ! proc
1110 105 : NULLIFY (real_msg)
1111 105 : NULLIFY (real_msg_gather)
1112 420 : msglen = SIZE(helium_env(1)%helium%pos(:, :, 1:helium_env(1)%helium%beads))
1113 315 : ALLOCATE (real_msg(msglen*helium_env(1)%helium%num_env))
1114 315 : ALLOCATE (real_msg_gather(msglen*helium_env(1)%helium%num_env))
1115 340329 : real_msg(:) = 0.0_dp
1116 240 : DO i = 1, SIZE(helium_env)
1117 172752 : real_msg((offset+i-1)*msglen+1:(offset+i)*msglen) = PACK(helium_env(i)%helium%pos(:, :, 1:helium_env(i)%helium%beads), .TRUE.)
1118 : END DO
1119 :
1120 : ! pass the message from all processors to logger%para_env%source
1121 680553 : CALL helium_env(1)%comm%sum(real_msg)
1122 680658 : real_msg_gather(:) = real_msg(:)
1123 :
1124 : ! update coordinates in the global input structure, only in
1125 : ! helium_env(1)
1126 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1127 : "MOTION%PINT%HELIUM%COORD%_DEFAULT_KEYWORD_", &
1128 105 : r_vals_ptr=real_msg_gather)
1129 :
1130 : ! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
1131 : ! assigned in section_vals_val_set - this memory will be used later on!
1132 : ! "The val becomes the owner of the array" - from section_vals_val_set docu
1133 105 : NULLIFY (real_msg_gather)
1134 :
1135 : ! DEALLOCATE since this array is only used locally
1136 105 : DEALLOCATE (real_msg)
1137 :
1138 : !
1139 : ! save permutation state
1140 : !
1141 : ! allocate the buffer for message passing
1142 105 : NULLIFY (int_msg_gather)
1143 105 : msglen = SIZE(helium_env(1)%helium%permutation)
1144 315 : ALLOCATE (int_msg_gather(msglen*helium_env(1)%helium%num_env))
1145 :
1146 : ! pass the message from all processors to logger%para_env%source
1147 5825 : int_msg_gather(:) = 0
1148 240 : DO i = 1, SIZE(helium_env)
1149 3150 : int_msg_gather((offset + i - 1)*msglen + 1:(offset + i)*msglen) = helium_env(i)%helium%permutation
1150 : END DO
1151 :
1152 11545 : CALL helium_env(1)%comm%sum(int_msg_gather)
1153 :
1154 : ! update permutation state in the global input structure
1155 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1156 : "MOTION%PINT%HELIUM%PERM%_DEFAULT_KEYWORD_", &
1157 105 : i_vals_ptr=int_msg_gather)
1158 :
1159 : ! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
1160 : ! assigned in section_vals_val_set - this memory will be used later on!
1161 : ! "The val becomes the owner of the array" - from section_vals_val_set docu
1162 105 : NULLIFY (int_msg_gather)
1163 :
1164 : !
1165 : ! save averages
1166 : !
1167 : ! update the weighting factor
1168 105 : itmp = helium_env(1)%helium%averages_iweight
1169 105 : IF (itmp < 0) THEN
1170 0 : itmp = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1171 : ELSE
1172 105 : itmp = itmp + helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1173 : END IF
1174 240 : DO i = 1, SIZE(helium_env)
1175 : CALL section_vals_val_set(helium_env(i)%helium%input, &
1176 : "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
1177 240 : i_val=itmp)
1178 : END DO
1179 :
1180 : ! allocate the buffer for message passing
1181 105 : NULLIFY (real_msg_gather)
1182 105 : msglen = 3
1183 315 : ALLOCATE (real_msg_gather(msglen*helium_env(1)%helium%num_env))
1184 :
1185 900 : real_msg_gather(:) = 0.0_dp
1186 : ! gather projected area from all processors
1187 240 : DO i = 1, SIZE(helium_env)
1188 645 : real_msg_gather((i - 1 + offset)*msglen + 1:(i + offset)*msglen) = helium_env(i)%helium%proarea%ravr(:)
1189 : END DO
1190 1695 : CALL helium_env(1)%comm%sum(real_msg_gather)
1191 :
1192 : ! update it in the global input structure
1193 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1194 : "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
1195 105 : r_vals_ptr=real_msg_gather)
1196 :
1197 : ! allocate the buffer for message passing
1198 105 : NULLIFY (real_msg_gather)
1199 105 : msglen = 3
1200 315 : ALLOCATE (real_msg_gather(msglen*helium_env(1)%helium%num_env))
1201 :
1202 900 : real_msg_gather(:) = 0.0_dp
1203 : ! gather projected area squared from all processors
1204 240 : DO i = 1, SIZE(helium_env)
1205 645 : real_msg_gather((i - 1 + offset)*msglen + 1:(i + offset)*msglen) = helium_env(i)%helium%prarea2%ravr(:)
1206 : END DO
1207 1695 : CALL helium_env(1)%comm%sum(real_msg_gather)
1208 :
1209 : ! update it in the global input structure
1210 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1211 : "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
1212 105 : r_vals_ptr=real_msg_gather)
1213 :
1214 : ! allocate the buffer for message passing
1215 105 : NULLIFY (real_msg_gather)
1216 105 : msglen = 3
1217 315 : ALLOCATE (real_msg_gather(msglen*helium_env(1)%helium%num_env))
1218 :
1219 900 : real_msg_gather(:) = 0.0_dp
1220 : ! gather winding number squared from all processors
1221 240 : DO i = 1, SIZE(helium_env)
1222 645 : real_msg_gather((i - 1 + offset)*msglen + 1:(i + offset)*msglen) = helium_env(i)%helium%wnmber2%ravr(:)
1223 : END DO
1224 1695 : CALL helium_env(1)%comm%sum(real_msg_gather)
1225 :
1226 : ! update it in the global input structure
1227 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1228 : "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
1229 105 : r_vals_ptr=real_msg_gather)
1230 :
1231 : ! allocate the buffer for message passing
1232 105 : NULLIFY (real_msg_gather)
1233 105 : msglen = 3
1234 315 : ALLOCATE (real_msg_gather(msglen*helium_env(1)%helium%num_env))
1235 :
1236 900 : real_msg_gather(:) = 0.0_dp
1237 : ! gather moment of inertia from all processors
1238 240 : DO i = 1, SIZE(helium_env)
1239 645 : real_msg_gather((i - 1 + offset)*msglen + 1:(i + offset)*msglen) = helium_env(i)%helium%mominer%ravr(:)
1240 : END DO
1241 1695 : CALL helium_env(1)%comm%sum(real_msg_gather)
1242 :
1243 : ! update it in the global input structure
1244 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1245 : "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
1246 105 : r_vals_ptr=real_msg_gather)
1247 :
1248 : ! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
1249 : ! assigned in section_vals_val_set - this memory will be used later on!
1250 : ! "The val becomes the owner of the array" - from section_vals_val_set docu
1251 105 : NULLIFY (real_msg_gather)
1252 :
1253 : !
1254 : ! save RNG state
1255 : !
1256 : ! pack RNG state on each processor to the local array and save in
1257 : ! gather with offset determined earlier
1258 : NULLIFY (real_msg)
1259 105 : msglen = 40
1260 105 : ALLOCATE (real_msg(msglen))
1261 : NULLIFY (real_msg_gather)
1262 315 : ALLOCATE (real_msg_gather(msglen*helium_env(1)%helium%num_env))
1263 10705 : real_msg_gather(:) = 0.0_dp
1264 :
1265 240 : DO i = 1, SIZE(helium_env)
1266 : CALL helium_env(i)%helium%rng_stream_uniform%get(bg=bg, cg=cg, ig=ig, &
1267 135 : buffer=bu, buffer_filled=lbf)
1268 135 : off = 0
1269 135 : real_msg(off + 1:off + 6) = PACK(bg, .TRUE.)
1270 135 : real_msg(off + 7:off + 12) = PACK(cg, .TRUE.)
1271 135 : real_msg(off + 13:off + 18) = PACK(ig, .TRUE.)
1272 135 : IF (lbf) THEN
1273 : bf = 1.0_dp
1274 : ELSE
1275 135 : bf = -1.0_dp
1276 : END IF
1277 135 : real_msg(off + 19) = bf
1278 135 : real_msg(off + 20) = bu
1279 : CALL helium_env(i)%helium%rng_stream_gaussian%get(bg=bg, cg=cg, ig=ig, &
1280 135 : buffer=bu, buffer_filled=lbf)
1281 135 : off = 20
1282 135 : real_msg(off + 1:off + 6) = PACK(bg, .TRUE.)
1283 135 : real_msg(off + 7:off + 12) = PACK(cg, .TRUE.)
1284 135 : real_msg(off + 13:off + 18) = PACK(ig, .TRUE.)
1285 135 : IF (lbf) THEN
1286 : bf = 1.0_dp
1287 : ELSE
1288 73 : bf = -1.0_dp
1289 : END IF
1290 135 : real_msg(off + 19) = bf
1291 135 : real_msg(off + 20) = bu
1292 :
1293 11175 : real_msg_gather((offset + i - 1)*msglen + 1:(offset + i)*msglen) = real_msg(:)
1294 : END DO
1295 :
1296 : ! Gather RNG state (in real_msg_gather vector) from all processors at
1297 : ! logger%para_env%source
1298 21305 : CALL helium_env(1)%comm%sum(real_msg_gather)
1299 :
1300 : ! update the RNG state in the global input structure
1301 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1302 : "MOTION%PINT%HELIUM%RNG_STATE%_DEFAULT_KEYWORD_", &
1303 105 : r_vals_ptr=real_msg_gather)
1304 :
1305 : ! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
1306 : ! assigned in section_vals_val_set - this memeory will be used later on!
1307 : ! "The val becomes the owner of the array" - from section_vals_val_set docu
1308 105 : NULLIFY (real_msg_gather)
1309 :
1310 : ! DEALLOCATE since this array is only used locally
1311 105 : DEALLOCATE (real_msg)
1312 :
1313 105 : IF (helium_env(1)%helium%solute_present) THEN
1314 : !
1315 : ! save forces on the solute
1316 : !
1317 : ! check that the number of values match the current runtime
1318 67 : reqlen = helium_env(1)%helium%solute_atoms*helium_env(1)%helium%solute_beads*3
1319 201 : msglen = SIZE(helium_env(1)%helium%force_avrg)
1320 67 : err_str = "Invalid size of HELIUM%FORCE: received '"
1321 67 : stmp = ""
1322 67 : WRITE (stmp, *) msglen
1323 : err_str = TRIM(ADJUSTL(err_str))// &
1324 67 : TRIM(ADJUSTL(stmp))//"' but expected '"
1325 67 : stmp = ""
1326 67 : WRITE (stmp, *) reqlen
1327 : err_str = TRIM(ADJUSTL(err_str))// &
1328 67 : TRIM(ADJUSTL(stmp))//"'."
1329 67 : IF (msgLEN /= reqlen) THEN
1330 0 : CPABORT(err_str)
1331 : END IF
1332 :
1333 : ! allocate the buffer to be saved and fill it with forces
1334 : ! forces should be the same on all processors, but we don't check that here
1335 67 : NULLIFY (real_msg_gather)
1336 201 : ALLOCATE (real_msg_gather(msglen))
1337 4531 : real_msg_gather(:) = PACK(helium_env(1)%helium%force_avrg, .TRUE.)
1338 :
1339 : ! update forces in the global input structure
1340 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1341 : "MOTION%PINT%HELIUM%FORCE%_DEFAULT_KEYWORD_", &
1342 67 : r_vals_ptr=real_msg_gather)
1343 :
1344 : ! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
1345 : ! assigned in section_vals_val_set - this memeory will be used later on!
1346 : ! "The val becomes the owner of the array" - from section_vals_val_set docu
1347 67 : NULLIFY (real_msg_gather)
1348 : END IF
1349 :
1350 : !
1351 : ! save the RDFs
1352 : !
1353 105 : IF (helium_env(1)%helium%rdf_present) THEN
1354 :
1355 : ! work on the temporary array so that accumulated data remains intact
1356 5010 : helium_env(1)%helium%rdf_inst(:, :) = 0.0_dp
1357 20 : DO i = 1, SIZE(helium_env)
1358 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :) + &
1359 5020 : helium_env(i)%helium%rdf_accu(:, :)
1360 : END DO
1361 :
1362 : ! average over processors / helium environments
1363 10010 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rdf_inst)
1364 10 : itmp = helium_env(1)%helium%num_env
1365 10 : invproc = 1.0_dp/REAL(itmp, dp)
1366 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*invproc
1367 :
1368 10 : nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1369 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps, dp)
1370 10 : iweight = helium_env(1)%helium%rdf_iweight
1371 : ! average over the old and the current density (observe the weights!)
1372 : helium_env(1)%helium%rdf_inst(:, :) = nsteps*helium_env(1)%helium%rdf_inst(:, :) + &
1373 5010 : iweight*helium_env(1)%helium%rdf_rstr(:, :)
1374 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps + iweight, dp)
1375 : ! update in the global input structure
1376 10 : NULLIFY (real_msg)
1377 30 : msglen = SIZE(helium_env(1)%helium%rdf_inst)
1378 30 : ALLOCATE (real_msg(msglen))
1379 2510 : real_msg(:) = PACK(helium_env(1)%helium%rdf_inst, .TRUE.)
1380 : CALL section_vals_val_set( &
1381 : helium_env(1)%helium%input, &
1382 : "MOTION%PINT%HELIUM%AVERAGES%RDF", &
1383 10 : r_vals_ptr=real_msg)
1384 10 : NULLIFY (real_msg)
1385 :
1386 : END IF
1387 :
1388 : !
1389 : ! save the densities
1390 : !
1391 105 : IF (helium_env(1)%helium%rho_present) THEN
1392 :
1393 : ! work on the temporary array so that accumulated data remains intact
1394 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = 0.0_dp
1395 20 : DO i = 1, SIZE(helium_env)
1396 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :) + &
1397 21120 : helium_env(i)%helium%rho_accu(:, :, :, :)
1398 : END DO
1399 :
1400 : ! average over processors / helium environments
1401 42210 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rho_inst)
1402 10 : itmp = helium_env(1)%helium%num_env
1403 10 : invproc = 1.0_dp/REAL(itmp, dp)
1404 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*invproc
1405 :
1406 10 : nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1407 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/REAL(nsteps, dp)
1408 10 : iweight = helium_env(1)%helium%averages_iweight
1409 : ! average over the old and the current density (observe the weights!)
1410 : helium_env(1)%helium%rho_inst(:, :, :, :) = nsteps*helium_env(1)%helium%rho_inst(:, :, :, :) + &
1411 21110 : iweight*helium_env(1)%helium%rho_rstr(:, :, :, :)
1412 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/REAL(nsteps + iweight, dp)
1413 :
1414 : ! update the densities in the global input structure
1415 10 : NULLIFY (real_msg)
1416 50 : msglen = SIZE(helium_env(1)%helium%rho_inst)
1417 30 : ALLOCATE (real_msg(msglen))
1418 10010 : real_msg(:) = PACK(helium_env(1)%helium%rho_inst, .TRUE.)
1419 : CALL section_vals_val_set( &
1420 : helium_env(1)%helium%input, &
1421 : "MOTION%PINT%HELIUM%AVERAGES%RHO", &
1422 10 : r_vals_ptr=real_msg)
1423 10 : NULLIFY (real_msg)
1424 :
1425 : END IF
1426 :
1427 : END IF ! ASSOCIATED(helium_env)
1428 :
1429 110 : CALL timestop(handle)
1430 :
1431 110 : END SUBROUTINE update_motion_helium
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief routine to dump thermostat CSVR energies
1435 : !> \param thermostat_energy ...
1436 : !> \param nsize ...
1437 : !> \param work_section ...
1438 : !> \par History
1439 : !> 10.2007 created [teo]
1440 : !> \author Teodoro Laino - University of Zurich
1441 : ! **************************************************************************************************
1442 422 : SUBROUTINE dump_csvr_energy_info(thermostat_energy, nsize, work_section)
1443 :
1444 : REAL(KIND=dp), DIMENSION(:), POINTER :: thermostat_energy
1445 : INTEGER, INTENT(IN) :: nsize
1446 : TYPE(section_vals_type), POINTER :: work_section
1447 :
1448 : INTEGER :: ik, irk, Nlist
1449 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
1450 : TYPE(section_type), POINTER :: section
1451 : TYPE(val_type), POINTER :: my_val, old_val
1452 :
1453 422 : CPASSERT(ASSOCIATED(work_section))
1454 422 : CPASSERT(work_section%ref_count > 0)
1455 :
1456 422 : NULLIFY (my_val, old_val, section, vals)
1457 :
1458 422 : section => work_section%section
1459 :
1460 422 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
1461 :
1462 422 : IF (ik == -2) THEN
1463 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
1464 0 : "_DEFAULT_KEYWORD_")
1465 : END IF
1466 :
1467 136 : DO
1468 558 : IF (SIZE(work_section%values, 2) == 1) EXIT
1469 136 : CALL section_vals_add_values(work_section)
1470 : END DO
1471 :
1472 422 : vals => work_section%values(ik, 1)%list
1473 422 : Nlist = 0
1474 :
1475 422 : IF (ASSOCIATED(vals)) THEN
1476 286 : Nlist = cp_sll_val_get_length(vals)
1477 : END IF
1478 :
1479 22798 : DO irk = 1, nsize
1480 22376 : CALL val_create(val=my_val, r_val=thermostat_energy(irk))
1481 22376 : IF (Nlist /= 0) THEN
1482 19024 : IF (irk == 1) THEN
1483 286 : new_pos => vals
1484 : ELSE
1485 18738 : new_pos => new_pos%rest
1486 : END IF
1487 19024 : old_val => new_pos%first_el
1488 19024 : CALL val_release(old_val)
1489 19024 : new_pos%first_el => my_val
1490 : ELSE
1491 3352 : IF (irk == 1) THEN
1492 136 : NULLIFY (new_pos)
1493 136 : CALL cp_sll_val_create(new_pos, first_el=my_val)
1494 136 : vals => new_pos
1495 : ELSE
1496 3216 : NULLIFY (new_pos%rest)
1497 3216 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
1498 3216 : new_pos => new_pos%rest
1499 : END IF
1500 : END IF
1501 22798 : NULLIFY (my_val)
1502 : END DO
1503 422 : work_section%values(ik, 1)%list => vals
1504 :
1505 422 : END SUBROUTINE dump_csvr_energy_info
1506 :
1507 : ! **************************************************************************************************
1508 : !> \brief Collect all information needed to dump the restart for CSVR
1509 : !> thermostat
1510 : !> \param csvr ...
1511 : !> \param para_env ...
1512 : !> \param csvr_section ...
1513 : !> \par History
1514 : !> 10.2007 created [tlaino] - University of Zurich
1515 : !> \author Teodoro Laino
1516 : ! **************************************************************************************************
1517 402 : SUBROUTINE dump_csvr_restart_info(csvr, para_env, csvr_section)
1518 :
1519 : TYPE(csvr_system_type), POINTER :: csvr
1520 : TYPE(mp_para_env_type), POINTER :: para_env
1521 : TYPE(section_vals_type), POINTER :: csvr_section
1522 :
1523 : CHARACTER(LEN=rng_record_length) :: rng_record
1524 : INTEGER :: i, my_index
1525 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dwork
1526 : REAL(KIND=dp) :: dum
1527 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: thermo_energy
1528 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
1529 : TYPE(section_vals_type), POINTER :: work_section
1530 :
1531 : ! Thermostat Energies
1532 :
1533 1206 : ALLOCATE (work(csvr%glob_num_csvr))
1534 :
1535 1206 : ALLOCATE (thermo_energy(csvr%loc_num_csvr))
1536 8504 : DO i = 1, csvr%loc_num_csvr
1537 8504 : thermo_energy(i) = csvr%nvt(i)%thermostat_energy
1538 : END DO
1539 : CALL get_kin_energies(csvr%map_info, csvr%loc_num_csvr, &
1540 : csvr%glob_num_csvr, thermo_energy, &
1541 402 : dum, para_env, array_kin=work)
1542 402 : DEALLOCATE (thermo_energy)
1543 :
1544 : ! If check passes then let's dump the info on the restart file
1545 402 : work_section => section_vals_get_subs_vals(csvr_section, "THERMOSTAT_ENERGY")
1546 402 : CALL dump_csvr_energy_info(work, csvr%glob_num_csvr, work_section)
1547 402 : DEALLOCATE (work)
1548 :
1549 : ! Thermostat Random Number info for restart
1550 402 : work_section => section_vals_get_subs_vals(csvr_section, "RNG_INIT")
1551 1206 : ALLOCATE (dwork(rng_record_length, csvr%glob_num_csvr))
1552 402 : dwork = 0
1553 8504 : DO i = 1, csvr%loc_num_csvr
1554 8102 : my_index = csvr%map_info%index(i)
1555 8102 : CALL csvr%nvt(i)%gaussian_rng_stream%dump(rng_record)
1556 8504 : CALL string_to_ascii(rng_record, dwork(:, my_index))
1557 : END DO
1558 :
1559 : ! Collect data if there was no communication in this thermostat
1560 402 : IF (csvr%map_info%dis_type == do_thermo_no_communication) THEN
1561 : ! Collect data if there was no communication in this thermostat
1562 148 : CALL para_env%sum(dwork)
1563 : ELSE
1564 : ! Perform some check and collect data in case of communicating thermostats
1565 254 : CALL communication_thermo_low2(dwork, rng_record_length, csvr%glob_num_csvr, para_env)
1566 : END IF
1567 402 : CALL section_rng_val_set(rng_section=work_section, nsize=csvr%glob_num_csvr, ascii=dwork)
1568 402 : DEALLOCATE (dwork)
1569 :
1570 804 : END SUBROUTINE dump_csvr_restart_info
1571 :
1572 : ! **************************************************************************************************
1573 : !> \brief Collect all information needed to dump the restart for AD_LANGEVIN
1574 : !> thermostat
1575 : !> \param al ...
1576 : !> \param para_env ...
1577 : !> \param al_section ...
1578 : !> \par History
1579 : !> 10.2007 created [tlaino] - University of Zurich
1580 : !> \author Teodoro Laino
1581 : ! **************************************************************************************************
1582 0 : SUBROUTINE dump_al_restart_info(al, para_env, al_section)
1583 :
1584 : TYPE(al_system_type), POINTER :: al
1585 : TYPE(mp_para_env_type), POINTER :: para_env
1586 : TYPE(section_vals_type), POINTER :: al_section
1587 :
1588 : INTEGER :: i
1589 : REAL(KIND=dp) :: dum
1590 : REAL(KIND=dp), DIMENSION(:), POINTER :: t_array, work
1591 : TYPE(section_vals_type), POINTER :: work_section
1592 :
1593 : ! chi and mass
1594 :
1595 0 : ALLOCATE (work(al%glob_num_al))
1596 0 : ALLOCATE (t_array(al%loc_num_al))
1597 :
1598 : ! copy chi into temporary
1599 0 : DO i = 1, al%loc_num_al
1600 0 : t_array(i) = al%nvt(i)%chi
1601 : END DO
1602 : ! consolidate into work
1603 : CALL get_kin_energies(al%map_info, al%loc_num_al, &
1604 : al%glob_num_al, t_array, &
1605 0 : dum, para_env, array_kin=work)
1606 :
1607 : ! If check passes then let's dump the info on the restart file
1608 0 : work_section => section_vals_get_subs_vals(al_section, "CHI")
1609 0 : CALL dump_csvr_energy_info(work, al%glob_num_al, work_section)
1610 :
1611 : ! copy mass into temporary
1612 0 : DO i = 1, al%loc_num_al
1613 0 : t_array(i) = al%nvt(i)%mass
1614 : END DO
1615 : ! consolidate into work
1616 : CALL get_kin_energies(al%map_info, al%loc_num_al, &
1617 : al%glob_num_al, t_array, &
1618 0 : dum, para_env, array_kin=work)
1619 :
1620 : ! If check passes then let's dump the info on the restart file
1621 0 : work_section => section_vals_get_subs_vals(al_section, "MASS")
1622 0 : CALL dump_csvr_energy_info(work, al%glob_num_al, work_section)
1623 :
1624 0 : DEALLOCATE (t_array)
1625 0 : DEALLOCATE (work)
1626 :
1627 0 : END SUBROUTINE dump_al_restart_info
1628 :
1629 : ! **************************************************************************************************
1630 : !> \brief Collect all information needed to dump the restart for GLE
1631 : !> thermostat
1632 : !> \param gle ...
1633 : !> \param para_env ...
1634 : !> \param gle_section ...
1635 : !> \author MI
1636 : ! **************************************************************************************************
1637 20 : SUBROUTINE dump_gle_restart_info(gle, para_env, gle_section)
1638 :
1639 : TYPE(gle_type), POINTER :: gle
1640 : TYPE(mp_para_env_type), POINTER :: para_env
1641 : TYPE(section_vals_type), POINTER :: gle_section
1642 :
1643 : CHARACTER(LEN=rng_record_length) :: rng_record
1644 : INTEGER :: counter, glob_num, i, iproc, j, loc_num
1645 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dwork
1646 20 : INTEGER, DIMENSION(:), POINTER :: gle_per_proc, index
1647 : REAL(dp) :: dum
1648 20 : REAL(dp), DIMENSION(:), POINTER :: s_tmp
1649 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: thermo_energy
1650 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
1651 : TYPE(section_vals_type), POINTER :: work_section
1652 :
1653 : ! Thermostat Energies
1654 :
1655 60 : ALLOCATE (work(gle%glob_num_gle))
1656 60 : ALLOCATE (thermo_energy(gle%loc_num_gle))
1657 3260 : DO i = 1, gle%loc_num_gle
1658 3260 : thermo_energy(i) = gle%nvt(i)%thermostat_energy
1659 : END DO
1660 : CALL get_kin_energies(gle%map_info, gle%loc_num_gle, &
1661 : gle%glob_num_gle, thermo_energy, &
1662 20 : dum, para_env, array_kin=work)
1663 20 : DEALLOCATE (thermo_energy)
1664 :
1665 : ! If check passes then let's dump the info on the restart file
1666 20 : work_section => section_vals_get_subs_vals(gle_section, "THERMOSTAT_ENERGY")
1667 20 : CALL dump_csvr_energy_info(work, gle%glob_num_gle, work_section)
1668 20 : DEALLOCATE (work)
1669 :
1670 : ! Thermostat Random Number info for restart
1671 20 : work_section => section_vals_get_subs_vals(gle_section, "RNG_INIT")
1672 20 : glob_num = gle%glob_num_gle
1673 20 : loc_num = gle%loc_num_gle
1674 60 : ALLOCATE (dwork(rng_record_length, glob_num))
1675 20 : dwork = 0
1676 3260 : DO i = 1, loc_num
1677 3240 : j = gle%map_info%index(i)
1678 3240 : CALL gle%nvt(i)%gaussian_rng_stream%dump(rng_record)
1679 3260 : CALL string_to_ascii(rng_record, dwork(:, j))
1680 : END DO
1681 :
1682 : ! Collect data if there was no communication in this thermostat
1683 20 : IF (gle%map_info%dis_type == do_thermo_no_communication) THEN
1684 : ! Collect data if there was no communication in this thermostat
1685 20 : CALL para_env%sum(dwork)
1686 : ELSE
1687 : ! Perform some check and collect data in case of communicating thermostats
1688 0 : CALL communication_thermo_low2(dwork, rng_record_length, glob_num, para_env)
1689 : END IF
1690 20 : CALL section_rng_val_set(rng_section=work_section, nsize=glob_num, ascii=dwork)
1691 20 : DEALLOCATE (dwork)
1692 :
1693 60 : ALLOCATE (gle_per_proc(para_env%num_pe))
1694 60 : gle_per_proc(:) = 0
1695 60 : CALL para_env%allgather(gle%loc_num_gle, gle_per_proc)
1696 :
1697 : ! Thermostat S variable info for restart
1698 20 : NULLIFY (s_tmp)
1699 60 : ALLOCATE (s_tmp((gle%ndim)*gle%glob_num_gle))
1700 32420 : s_tmp = 0.0_dp
1701 :
1702 20 : NULLIFY (work, index)
1703 60 : DO iproc = 1, para_env%num_pe
1704 40 : CALL reallocate(work, 1, gle_per_proc(iproc)*(gle%ndim))
1705 40 : CALL reallocate(index, 1, gle_per_proc(iproc))
1706 40 : IF (para_env%mepos == (iproc - 1)) THEN
1707 3260 : INDEX(:) = 0
1708 20 : counter = 0
1709 120 : DO i = 1, gle%ndim
1710 16320 : DO j = 1, gle%loc_num_gle
1711 16200 : counter = counter + 1
1712 16200 : work(counter) = gle%nvt(j)%s(i)
1713 16300 : INDEX(j) = gle%map_info%index(j)
1714 : END DO
1715 : END DO
1716 : ELSE
1717 16220 : work(:) = 0.0_dp
1718 : END IF
1719 64840 : CALL para_env%bcast(work, iproc - 1)
1720 13000 : CALL para_env%bcast(index, iproc - 1)
1721 40 : counter = 0
1722 260 : DO i = 1, gle%ndim
1723 32640 : DO j = 1, gle_per_proc(iproc)
1724 32400 : counter = counter + 1
1725 32600 : s_tmp((INDEX(j) - 1)*(gle%ndim) + i) = work(counter)
1726 : END DO
1727 : END DO
1728 : END DO
1729 :
1730 20 : IF (SIZE(s_tmp) > 0) THEN
1731 20 : work_section => section_vals_get_subs_vals(gle_section, "S")
1732 20 : CALL section_vals_val_set(work_section, "_DEFAULT_KEYWORD_", r_vals_ptr=s_tmp)
1733 : ELSE
1734 0 : DEALLOCATE (s_tmp)
1735 : END IF
1736 :
1737 20 : DEALLOCATE (gle_per_proc)
1738 20 : DEALLOCATE (work)
1739 20 : DEALLOCATE (index)
1740 :
1741 40 : END SUBROUTINE dump_gle_restart_info
1742 :
1743 : ! **************************************************************************************************
1744 : !> \brief Collect all information needed to dump the restart for Nose-Hoover
1745 : !> thermostat
1746 : !> \param nhc ...
1747 : !> \param para_env ...
1748 : !> \param eta ...
1749 : !> \param veta ...
1750 : !> \param fnhc ...
1751 : !> \param mnhc ...
1752 : !> \par History
1753 : !> 10.2007 created [tlaino] - University of Zurich
1754 : !> \author Teodoro Laino
1755 : ! **************************************************************************************************
1756 964 : SUBROUTINE collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
1757 :
1758 : TYPE(lnhc_parameters_type), POINTER :: nhc
1759 : TYPE(mp_para_env_type), POINTER :: para_env
1760 : REAL(KIND=dp), DIMENSION(:), POINTER :: eta, veta, fnhc, mnhc
1761 :
1762 : INTEGER :: counter, i, iproc, j, nhc_len, num_nhc, &
1763 : numneed, tot_nhcneed
1764 964 : INTEGER, DIMENSION(:), POINTER :: index, nhc_per_proc
1765 964 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
1766 : TYPE(map_info_type), POINTER :: map_info
1767 :
1768 964 : nhc_len = SIZE(nhc%nvt, 1)
1769 964 : num_nhc = nhc%loc_num_nhc
1770 964 : numneed = num_nhc
1771 964 : map_info => nhc%map_info
1772 2892 : ALLOCATE (nhc_per_proc(para_env%num_pe))
1773 2892 : nhc_per_proc(:) = 0
1774 :
1775 2892 : CALL para_env%allgather(numneed, nhc_per_proc)
1776 964 : tot_nhcneed = nhc%glob_num_nhc
1777 :
1778 964 : NULLIFY (work, index)
1779 : !-----------------------------------------------------------------------------
1780 : !-----------------------------------------------------------------------------
1781 : ! nhc%eta
1782 : !-----------------------------------------------------------------------------
1783 2892 : ALLOCATE (eta(tot_nhcneed*nhc_len))
1784 2892 : DO iproc = 1, para_env%num_pe
1785 1928 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1786 1928 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1787 1928 : IF (para_env%mepos == (iproc - 1)) THEN
1788 68943 : INDEX(:) = 0
1789 : counter = 0
1790 4670 : DO i = 1, nhc_len
1791 257195 : DO j = 1, num_nhc
1792 252525 : counter = counter + 1
1793 252525 : work(counter) = nhc%nvt(i, j)%eta
1794 256231 : INDEX(j) = map_info%index(j)
1795 : END DO
1796 : END DO
1797 : ELSE
1798 253489 : work(:) = 0.0_dp
1799 : END IF
1800 1012028 : CALL para_env%bcast(work, iproc - 1)
1801 273844 : CALL para_env%bcast(index, iproc - 1)
1802 1928 : counter = 0
1803 10304 : DO i = 1, nhc_len
1804 514390 : DO j = 1, nhc_per_proc(iproc)
1805 505050 : counter = counter + 1
1806 512462 : eta((INDEX(j) - 1)*nhc_len + i) = work(counter)
1807 : END DO
1808 : END DO
1809 : END DO
1810 : !-----------------------------------------------------------------------------
1811 : !-----------------------------------------------------------------------------
1812 : ! nhc%veta
1813 : !-----------------------------------------------------------------------------
1814 1928 : ALLOCATE (veta(tot_nhcneed*nhc_len))
1815 2892 : DO iproc = 1, para_env%num_pe
1816 1928 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1817 1928 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1818 1928 : IF (para_env%mepos == (iproc - 1)) THEN
1819 68943 : INDEX(:) = 0
1820 : counter = 0
1821 4670 : DO i = 1, nhc_len
1822 257195 : DO j = 1, num_nhc
1823 252525 : counter = counter + 1
1824 252525 : work(counter) = nhc%nvt(i, j)%v
1825 256231 : INDEX(j) = map_info%index(j)
1826 : END DO
1827 : END DO
1828 : ELSE
1829 253489 : work(:) = 0.0_dp
1830 : END IF
1831 1012028 : CALL para_env%bcast(work, iproc - 1)
1832 273844 : CALL para_env%bcast(index, iproc - 1)
1833 1928 : counter = 0
1834 10304 : DO i = 1, nhc_len
1835 514390 : DO j = 1, nhc_per_proc(iproc)
1836 505050 : counter = counter + 1
1837 512462 : veta((INDEX(j) - 1)*nhc_len + i) = work(counter)
1838 : END DO
1839 : END DO
1840 : END DO
1841 : !-----------------------------------------------------------------------------
1842 : !-----------------------------------------------------------------------------
1843 : ! nhc%force
1844 : !-----------------------------------------------------------------------------
1845 1928 : ALLOCATE (fnhc(tot_nhcneed*nhc_len))
1846 2892 : DO iproc = 1, para_env%num_pe
1847 1928 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1848 1928 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1849 1928 : IF (para_env%mepos == (iproc - 1)) THEN
1850 68943 : INDEX(:) = 0
1851 : counter = 0
1852 4670 : DO i = 1, nhc_len
1853 257195 : DO j = 1, num_nhc
1854 252525 : counter = counter + 1
1855 252525 : work(counter) = nhc%nvt(i, j)%f
1856 256231 : INDEX(j) = map_info%index(j)
1857 : END DO
1858 : END DO
1859 : ELSE
1860 253489 : work(:) = 0.0_dp
1861 : END IF
1862 1012028 : CALL para_env%bcast(work, iproc - 1)
1863 273844 : CALL para_env%bcast(index, iproc - 1)
1864 1928 : counter = 0
1865 10304 : DO i = 1, nhc_len
1866 514390 : DO j = 1, nhc_per_proc(iproc)
1867 505050 : counter = counter + 1
1868 512462 : fnhc((INDEX(j) - 1)*nhc_len + i) = work(counter)
1869 : END DO
1870 : END DO
1871 : END DO
1872 : !-----------------------------------------------------------------------------
1873 : !-----------------------------------------------------------------------------
1874 : ! nhc%mass
1875 : !-----------------------------------------------------------------------------
1876 1928 : ALLOCATE (mnhc(tot_nhcneed*nhc_len))
1877 2892 : DO iproc = 1, para_env%num_pe
1878 1928 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1879 1928 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1880 1928 : IF (para_env%mepos == (iproc - 1)) THEN
1881 68943 : INDEX(:) = 0
1882 : counter = 0
1883 4670 : DO i = 1, nhc_len
1884 257195 : DO j = 1, num_nhc
1885 252525 : counter = counter + 1
1886 252525 : work(counter) = nhc%nvt(i, j)%mass
1887 256231 : INDEX(j) = map_info%index(j)
1888 : END DO
1889 : END DO
1890 : ELSE
1891 253489 : work(:) = 0.0_dp
1892 : END IF
1893 1012028 : CALL para_env%bcast(work, iproc - 1)
1894 273844 : CALL para_env%bcast(index, iproc - 1)
1895 1928 : counter = 0
1896 10304 : DO i = 1, nhc_len
1897 514390 : DO j = 1, nhc_per_proc(iproc)
1898 505050 : counter = counter + 1
1899 512462 : mnhc((INDEX(j) - 1)*nhc_len + i) = work(counter)
1900 : END DO
1901 : END DO
1902 : END DO
1903 :
1904 964 : DEALLOCATE (work)
1905 964 : DEALLOCATE (index)
1906 964 : DEALLOCATE (nhc_per_proc)
1907 :
1908 964 : END SUBROUTINE collect_nose_restart_info
1909 :
1910 : ! **************************************************************************************************
1911 : !> \brief routine to dump NEB coordinates and velocities section.. fast implementation
1912 : !> \param coord_section ...
1913 : !> \param array ...
1914 : !> \param narray ...
1915 : !> \param nsize ...
1916 : !> \param nfield ...
1917 : !> \param particle_set ...
1918 : !> \param conv_factor ...
1919 : !> \par History
1920 : !> 12.2006 created [teo]
1921 : !> \author Teodoro Laino
1922 : ! **************************************************************************************************
1923 7148 : SUBROUTINE section_neb_coord_val_set(coord_section, array, narray, nsize, nfield, &
1924 : particle_set, conv_factor)
1925 :
1926 : TYPE(section_vals_type), POINTER :: coord_section
1927 : REAL(KIND=dp), DIMENSION(*) :: array
1928 : INTEGER, INTENT(IN) :: narray, nsize, nfield
1929 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1930 : REAL(KIND=dp) :: conv_factor
1931 :
1932 : INTEGER :: ik, irk, Nlist
1933 7148 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_c
1934 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
1935 : TYPE(section_type), POINTER :: section
1936 : TYPE(val_type), POINTER :: my_val, old_val
1937 :
1938 7148 : NULLIFY (my_val, old_val, section, vals)
1939 0 : CPASSERT(ASSOCIATED(coord_section))
1940 7148 : CPASSERT(coord_section%ref_count > 0)
1941 7148 : section => coord_section%section
1942 7148 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
1943 7148 : IF (ik == -2) THEN
1944 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
1945 0 : "_DEFAULT_KEYWORD_")
1946 : END IF
1947 364 : DO
1948 7512 : IF (SIZE(coord_section%values, 2) == 1) EXIT
1949 364 : CALL section_vals_add_values(coord_section)
1950 : END DO
1951 7148 : vals => coord_section%values(ik, 1)%list
1952 7148 : Nlist = 0
1953 7148 : IF (ASSOCIATED(vals)) THEN
1954 6784 : Nlist = cp_sll_val_get_length(vals)
1955 : END IF
1956 270192 : DO irk = 1, nsize/nfield
1957 789132 : ALLOCATE (my_c(nfield))
1958 263044 : IF (nfield == 3) THEN
1959 1051696 : my_c(1:3) = get_particle_pos_or_vel(irk, particle_set, array(1:narray))
1960 1051696 : my_c(1:3) = my_c(1:3)*conv_factor
1961 : ELSE
1962 120 : my_c(1) = array(irk)
1963 : END IF
1964 263044 : CALL val_create(my_val, r_vals_ptr=my_c)
1965 :
1966 263044 : IF (Nlist /= 0) THEN
1967 241140 : IF (irk == 1) THEN
1968 6784 : new_pos => vals
1969 : ELSE
1970 234356 : new_pos => new_pos%rest
1971 : END IF
1972 241140 : old_val => new_pos%first_el
1973 241140 : CALL val_release(old_val)
1974 241140 : new_pos%first_el => my_val
1975 : ELSE
1976 21904 : IF (irk == 1) THEN
1977 364 : NULLIFY (new_pos)
1978 364 : CALL cp_sll_val_create(new_pos, first_el=my_val)
1979 364 : vals => new_pos
1980 : ELSE
1981 21540 : NULLIFY (new_pos%rest)
1982 21540 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
1983 21540 : new_pos => new_pos%rest
1984 : END IF
1985 : END IF
1986 270192 : NULLIFY (my_val)
1987 : END DO
1988 :
1989 7148 : coord_section%values(ik, 1)%list => vals
1990 :
1991 7148 : END SUBROUTINE section_neb_coord_val_set
1992 :
1993 : ! **************************************************************************************************
1994 : !> \brief Set the nose structure like restart
1995 : !> \param work_section ...
1996 : !> \param eta ...
1997 : !> \param veta ...
1998 : !> \param fnhc ...
1999 : !> \param mnhc ...
2000 : !> \par History
2001 : !> 01.2006 created [teo]
2002 : !> \author Teodoro Laino
2003 : ! **************************************************************************************************
2004 1414 : SUBROUTINE set_template_restart(work_section, eta, veta, fnhc, mnhc)
2005 :
2006 : TYPE(section_vals_type), POINTER :: work_section
2007 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eta, veta, fnhc, mnhc
2008 :
2009 : TYPE(section_vals_type), POINTER :: coord, force, mass, velocity
2010 :
2011 1414 : NULLIFY (coord, force, velocity, mass)
2012 1414 : IF (PRESENT(eta)) THEN
2013 1016 : IF (SIZE(eta) > 0) THEN
2014 1016 : coord => section_vals_get_subs_vals(work_section, "COORD")
2015 1016 : CALL section_vals_val_set(coord, "_DEFAULT_KEYWORD_", r_vals_ptr=eta)
2016 : ELSE
2017 0 : DEALLOCATE (eta)
2018 : END IF
2019 : END IF
2020 1414 : IF (PRESENT(veta)) THEN
2021 1414 : IF (SIZE(veta) > 0) THEN
2022 1414 : velocity => section_vals_get_subs_vals(work_section, "VELOCITY")
2023 1414 : CALL section_vals_val_set(velocity, "_DEFAULT_KEYWORD_", r_vals_ptr=veta)
2024 : ELSE
2025 0 : DEALLOCATE (veta)
2026 : END IF
2027 : END IF
2028 1414 : IF (PRESENT(fnhc)) THEN
2029 1016 : IF (SIZE(fnhc) > 0) THEN
2030 1016 : force => section_vals_get_subs_vals(work_section, "FORCE")
2031 1016 : CALL section_vals_val_set(force, "_DEFAULT_KEYWORD_", r_vals_ptr=fnhc)
2032 : ELSE
2033 0 : DEALLOCATE (fnhc)
2034 : END IF
2035 : END IF
2036 1414 : IF (PRESENT(mnhc)) THEN
2037 1414 : IF (SIZE(mnhc) > 0) THEN
2038 1414 : mass => section_vals_get_subs_vals(work_section, "MASS")
2039 1414 : CALL section_vals_val_set(mass, "_DEFAULT_KEYWORD_", r_vals_ptr=mnhc)
2040 : ELSE
2041 0 : DEALLOCATE (mnhc)
2042 : END IF
2043 : END IF
2044 :
2045 1414 : END SUBROUTINE set_template_restart
2046 :
2047 : ! **************************************************************************************************
2048 : !> \brief routine to dump hills information during metadynamics run
2049 : !> \param ss_section ...
2050 : !> \param meta_env ...
2051 : !> \par History
2052 : !> 02.2006 created [teo]
2053 : !> \author Teodoro Laino
2054 : ! **************************************************************************************************
2055 784 : SUBROUTINE meta_hills_val_set_ss(ss_section, meta_env)
2056 :
2057 : TYPE(section_vals_type), POINTER :: ss_section
2058 : TYPE(meta_env_type), POINTER :: meta_env
2059 :
2060 : INTEGER :: ik, irk, lsize, Nlist
2061 784 : REAL(KIND=dp), DIMENSION(:), POINTER :: ss_val
2062 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2063 : TYPE(section_type), POINTER :: section
2064 : TYPE(val_type), POINTER :: my_val, old_val
2065 :
2066 784 : NULLIFY (my_val, old_val, section, vals)
2067 0 : CPASSERT(ASSOCIATED(ss_section))
2068 784 : CPASSERT(ss_section%ref_count > 0)
2069 784 : section => ss_section%section
2070 784 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2071 784 : IF (ik == -2) THEN
2072 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2073 0 : "_DEFAULT_KEYWORD_")
2074 : END IF
2075 100 : DO
2076 884 : IF (SIZE(ss_section%values, 2) == 1) EXIT
2077 100 : CALL section_vals_add_values(ss_section)
2078 : END DO
2079 784 : vals => ss_section%values(ik, 1)%list
2080 784 : Nlist = 0
2081 784 : IF (ASSOCIATED(vals)) THEN
2082 684 : Nlist = cp_sll_val_get_length(vals)
2083 : END IF
2084 784 : lsize = SIZE(meta_env%hills_env%ss_history, 1)
2085 12938 : DO irk = 1, meta_env%hills_env%n_hills
2086 36462 : ALLOCATE (ss_val(lsize))
2087 : ! Always stored in A.U.
2088 49256 : ss_val = meta_env%hills_env%ss_history(:, irk)
2089 12154 : CALL val_create(my_val, r_vals_ptr=ss_val)
2090 :
2091 12154 : IF (irk <= Nlist) THEN
2092 10980 : IF (irk == 1) THEN
2093 684 : new_pos => vals
2094 : ELSE
2095 10296 : new_pos => new_pos%rest
2096 : END IF
2097 10980 : old_val => new_pos%first_el
2098 10980 : CALL val_release(old_val)
2099 10980 : new_pos%first_el => my_val
2100 : ELSE
2101 1174 : IF (irk == 1) THEN
2102 100 : NULLIFY (new_pos)
2103 100 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2104 100 : vals => new_pos
2105 : ELSE
2106 1074 : NULLIFY (new_pos%rest)
2107 1074 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2108 1074 : new_pos => new_pos%rest
2109 : END IF
2110 : END IF
2111 12938 : NULLIFY (my_val)
2112 : END DO
2113 :
2114 784 : ss_section%values(ik, 1)%list => vals
2115 :
2116 784 : END SUBROUTINE meta_hills_val_set_ss
2117 :
2118 : ! **************************************************************************************************
2119 : !> \brief routine to dump hills information during metadynamics run
2120 : !> \param ds_section ...
2121 : !> \param meta_env ...
2122 : !> \par History
2123 : !> 02.2006 created [teo]
2124 : !> \author Teodoro Laino
2125 : ! **************************************************************************************************
2126 784 : SUBROUTINE meta_hills_val_set_ds(ds_section, meta_env)
2127 :
2128 : TYPE(section_vals_type), POINTER :: ds_section
2129 : TYPE(meta_env_type), POINTER :: meta_env
2130 :
2131 : INTEGER :: ik, irk, lsize, Nlist
2132 784 : REAL(KIND=dp), DIMENSION(:), POINTER :: ds_val
2133 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2134 : TYPE(section_type), POINTER :: section
2135 : TYPE(val_type), POINTER :: my_val, old_val
2136 :
2137 784 : NULLIFY (my_val, old_val, section, vals)
2138 0 : CPASSERT(ASSOCIATED(ds_section))
2139 784 : CPASSERT(ds_section%ref_count > 0)
2140 784 : section => ds_section%section
2141 784 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2142 784 : IF (ik == -2) THEN
2143 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2144 0 : "_DEFAULT_KEYWORD_")
2145 : END IF
2146 100 : DO
2147 884 : IF (SIZE(ds_section%values, 2) == 1) EXIT
2148 100 : CALL section_vals_add_values(ds_section)
2149 : END DO
2150 784 : vals => ds_section%values(ik, 1)%list
2151 784 : Nlist = 0
2152 784 : IF (ASSOCIATED(vals)) THEN
2153 684 : Nlist = cp_sll_val_get_length(vals)
2154 : END IF
2155 784 : lsize = SIZE(meta_env%hills_env%delta_s_history, 1)
2156 12938 : DO irk = 1, meta_env%hills_env%n_hills
2157 36462 : ALLOCATE (ds_val(lsize))
2158 : ! Always stored in A.U.
2159 49256 : ds_val = meta_env%hills_env%delta_s_history(:, irk)
2160 12154 : CALL val_create(my_val, r_vals_ptr=ds_val)
2161 :
2162 12154 : IF (irk <= Nlist) THEN
2163 10980 : IF (irk == 1) THEN
2164 684 : new_pos => vals
2165 : ELSE
2166 10296 : new_pos => new_pos%rest
2167 : END IF
2168 10980 : old_val => new_pos%first_el
2169 10980 : CALL val_release(old_val)
2170 10980 : new_pos%first_el => my_val
2171 : ELSE
2172 1174 : IF (irk == 1) THEN
2173 100 : NULLIFY (new_pos)
2174 100 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2175 100 : vals => new_pos
2176 : ELSE
2177 1074 : NULLIFY (new_pos%rest)
2178 1074 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2179 1074 : new_pos => new_pos%rest
2180 : END IF
2181 : END IF
2182 12938 : NULLIFY (my_val)
2183 : END DO
2184 :
2185 784 : ds_section%values(ik, 1)%list => vals
2186 :
2187 784 : END SUBROUTINE meta_hills_val_set_ds
2188 :
2189 : ! **************************************************************************************************
2190 : !> \brief routine to dump hills information during metadynamics run
2191 : !> \param ww_section ...
2192 : !> \param meta_env ...
2193 : !> \par History
2194 : !> 02.2006 created [teo]
2195 : !> \author Teodoro Laino
2196 : ! **************************************************************************************************
2197 784 : SUBROUTINE meta_hills_val_set_ww(ww_section, meta_env)
2198 :
2199 : TYPE(section_vals_type), POINTER :: ww_section
2200 : TYPE(meta_env_type), POINTER :: meta_env
2201 :
2202 : INTEGER :: ik, irk, lsize, Nlist
2203 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2204 : TYPE(section_type), POINTER :: section
2205 : TYPE(val_type), POINTER :: my_val, old_val
2206 :
2207 784 : NULLIFY (my_val, old_val, section, vals)
2208 784 : CPASSERT(ASSOCIATED(ww_section))
2209 784 : CPASSERT(ww_section%ref_count > 0)
2210 784 : section => ww_section%section
2211 784 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2212 784 : IF (ik == -2) THEN
2213 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2214 0 : "_DEFAULT_KEYWORD_")
2215 : END IF
2216 100 : DO
2217 884 : IF (SIZE(ww_section%values, 2) == 1) EXIT
2218 100 : CALL section_vals_add_values(ww_section)
2219 : END DO
2220 784 : vals => ww_section%values(ik, 1)%list
2221 784 : Nlist = 0
2222 784 : IF (ASSOCIATED(vals)) THEN
2223 684 : Nlist = cp_sll_val_get_length(vals)
2224 : END IF
2225 784 : lsize = meta_env%hills_env%n_hills
2226 12938 : DO irk = 1, lsize
2227 12154 : CALL val_create(my_val, r_val=meta_env%hills_env%ww_history(irk))
2228 :
2229 12154 : IF (irk <= Nlist) THEN
2230 10980 : IF (irk == 1) THEN
2231 684 : new_pos => vals
2232 : ELSE
2233 10296 : new_pos => new_pos%rest
2234 : END IF
2235 10980 : old_val => new_pos%first_el
2236 10980 : CALL val_release(old_val)
2237 10980 : new_pos%first_el => my_val
2238 : ELSE
2239 1174 : IF (irk == 1) THEN
2240 100 : NULLIFY (new_pos)
2241 100 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2242 100 : vals => new_pos
2243 : ELSE
2244 1074 : NULLIFY (new_pos%rest)
2245 1074 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2246 1074 : new_pos => new_pos%rest
2247 : END IF
2248 : END IF
2249 12938 : NULLIFY (my_val)
2250 : END DO
2251 :
2252 784 : ww_section%values(ik, 1)%list => vals
2253 :
2254 784 : END SUBROUTINE meta_hills_val_set_ww
2255 :
2256 : ! **************************************************************************************************
2257 : !> \brief routine to dump hills information during metadynamics run
2258 : !> \param invdt_section ...
2259 : !> \param meta_env ...
2260 : !> \par History
2261 : !> 12.2009 created [seb]
2262 : !> \author SC
2263 : ! **************************************************************************************************
2264 2 : SUBROUTINE meta_hills_val_set_dt(invdt_section, meta_env)
2265 :
2266 : TYPE(section_vals_type), POINTER :: invdt_section
2267 : TYPE(meta_env_type), POINTER :: meta_env
2268 :
2269 : INTEGER :: ik, irk, lsize, Nlist
2270 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2271 : TYPE(section_type), POINTER :: section
2272 : TYPE(val_type), POINTER :: my_val, old_val
2273 :
2274 2 : NULLIFY (my_val, old_val, section, vals)
2275 2 : CPASSERT(ASSOCIATED(invdt_section))
2276 2 : CPASSERT(invdt_section%ref_count > 0)
2277 2 : section => invdt_section%section
2278 2 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2279 2 : IF (ik == -2) THEN
2280 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2281 0 : "_DEFAULT_KEYWORD_")
2282 : END IF
2283 2 : DO
2284 4 : IF (SIZE(invdt_section%values, 2) == 1) EXIT
2285 2 : CALL section_vals_add_values(invdt_section)
2286 : END DO
2287 2 : vals => invdt_section%values(ik, 1)%list
2288 2 : Nlist = 0
2289 2 : IF (ASSOCIATED(vals)) THEN
2290 0 : Nlist = cp_sll_val_get_length(vals)
2291 : END IF
2292 2 : lsize = meta_env%hills_env%n_hills
2293 6 : DO irk = 1, lsize
2294 4 : CALL val_create(my_val, r_val=meta_env%hills_env%invdt_history(irk))
2295 :
2296 4 : IF (irk <= Nlist) THEN
2297 0 : IF (irk == 1) THEN
2298 0 : new_pos => vals
2299 : ELSE
2300 0 : new_pos => new_pos%rest
2301 : END IF
2302 0 : old_val => new_pos%first_el
2303 0 : CALL val_release(old_val)
2304 0 : new_pos%first_el => my_val
2305 : ELSE
2306 4 : IF (irk == 1) THEN
2307 2 : NULLIFY (new_pos)
2308 2 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2309 2 : vals => new_pos
2310 : ELSE
2311 2 : NULLIFY (new_pos%rest)
2312 2 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2313 2 : new_pos => new_pos%rest
2314 : END IF
2315 : END IF
2316 6 : NULLIFY (my_val)
2317 : END DO
2318 2 : invdt_section%values(ik, 1)%list => vals
2319 2 : END SUBROUTINE meta_hills_val_set_dt
2320 :
2321 : ! **************************************************************************************************
2322 : !> \brief Write all input sections scaling in size with the number of atoms
2323 : !> in the system to an external file in binary format
2324 : !> \param output_unit binary file to write to
2325 : !> \param log_unit unit for logging debug information
2326 : !> \param root_section ...
2327 : !> \param md_env ...
2328 : !> \param force_env ...
2329 : !> \par History
2330 : !> - Creation (10.02.2011,MK)
2331 : !> \author Matthias Krack (MK)
2332 : !> \version 1.0
2333 : ! **************************************************************************************************
2334 272 : SUBROUTINE write_binary_restart(output_unit, log_unit, root_section, md_env, force_env)
2335 :
2336 : INTEGER, INTENT(IN) :: output_unit, log_unit
2337 : TYPE(section_vals_type), POINTER :: root_section
2338 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
2339 : TYPE(force_env_type), OPTIONAL, POINTER :: force_env
2340 :
2341 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_binary_restart'
2342 :
2343 : CHARACTER(LEN=default_path_length) :: binary_restart_file_name
2344 : CHARACTER(LEN=default_string_length) :: section_label
2345 : INTEGER :: handle, iatom, icore, ikind, imolecule, ishell, istat, n_char_size, n_dp_size, &
2346 : n_int_size, natom, natomkind, ncore, nhc_size, nmolecule, nmoleculekind, nshell, &
2347 : print_level, run_type
2348 272 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf, imol
2349 : LOGICAL :: print_info, write_velocities
2350 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rbuf
2351 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2352 : TYPE(cp_subsys_type), POINTER :: subsys
2353 : TYPE(force_env_type), POINTER :: my_force_env
2354 : TYPE(lnhc_parameters_type), POINTER :: nhc
2355 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2356 : TYPE(molecule_list_type), POINTER :: molecules
2357 : TYPE(mp_para_env_type), POINTER :: para_env
2358 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
2359 : shell_particles
2360 : TYPE(thermostat_type), POINTER :: thermostat_part, thermostat_shell
2361 :
2362 272 : CALL timeset(routineN, handle)
2363 :
2364 272 : NULLIFY (atomic_kinds)
2365 272 : NULLIFY (core_particles)
2366 272 : NULLIFY (molecule_kinds)
2367 272 : NULLIFY (molecules)
2368 272 : NULLIFY (my_force_env)
2369 272 : NULLIFY (para_env)
2370 272 : NULLIFY (particles)
2371 272 : NULLIFY (shell_particles)
2372 272 : NULLIFY (subsys)
2373 272 : NULLIFY (thermostat_part)
2374 272 : NULLIFY (thermostat_shell)
2375 :
2376 272 : IF (PRESENT(md_env)) THEN
2377 : CALL get_md_env(md_env=md_env, &
2378 : force_env=my_force_env, &
2379 : thermostat_part=thermostat_part, &
2380 272 : thermostat_shell=thermostat_shell)
2381 0 : ELSE IF (PRESENT(force_env)) THEN
2382 0 : my_force_env => force_env
2383 : END IF
2384 :
2385 272 : IF (.NOT. ASSOCIATED(my_force_env)) THEN
2386 0 : CALL timestop(handle)
2387 0 : RETURN
2388 : END IF
2389 :
2390 272 : CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", i_val=print_level)
2391 :
2392 272 : IF (print_level > 1) THEN
2393 272 : print_info = .TRUE.
2394 : ELSE
2395 0 : print_info = .FALSE.
2396 : END IF
2397 :
2398 272 : CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=run_type)
2399 : write_velocities = ((run_type == mol_dyn_run) .OR. &
2400 : (run_type == mon_car_run) .OR. &
2401 272 : (run_type == pint_run))
2402 :
2403 : CALL force_env_get(force_env=my_force_env, &
2404 : para_env=para_env, &
2405 272 : subsys=subsys)
2406 : CALL cp_subsys_get(subsys, &
2407 : atomic_kinds=atomic_kinds, &
2408 : particles=particles, &
2409 : natom=natom, &
2410 : core_particles=core_particles, &
2411 : ncore=ncore, &
2412 : shell_particles=shell_particles, &
2413 : nshell=nshell, &
2414 : molecule_kinds=molecule_kinds, &
2415 272 : molecules=molecules)
2416 :
2417 272 : natomkind = atomic_kinds%n_els
2418 272 : IF (ASSOCIATED(molecule_kinds)) THEN
2419 272 : nmoleculekind = molecule_kinds%n_els
2420 : ELSE
2421 0 : nmoleculekind = 0
2422 : END IF
2423 :
2424 272 : IF (ASSOCIATED(molecules)) THEN
2425 272 : nmolecule = molecules%n_els
2426 : ELSE
2427 0 : nmolecule = 0
2428 : END IF
2429 :
2430 272 : n_char_size = 0 ! init
2431 272 : n_int_size = 0 ! init
2432 272 : n_dp_size = 0 ! init
2433 :
2434 272 : IF (output_unit > 0) THEN ! only ionode
2435 :
2436 136 : IF (print_info) THEN
2437 136 : INQUIRE (UNIT=output_unit, NAME=binary_restart_file_name, IOSTAT=istat)
2438 136 : IF (istat /= 0) THEN
2439 : CALL cp_abort(__LOCATION__, &
2440 : "An error occurred inquiring logical unit <"// &
2441 : TRIM(ADJUSTL(cp_to_string(output_unit)))// &
2442 0 : "> which should be linked to the binary restart file")
2443 : END IF
2444 136 : IF (log_unit > 0) THEN
2445 : WRITE (UNIT=log_unit, FMT="(T2,A,/,/,(T3,A,T71,I10))") &
2446 136 : "Writing binary restart file "//TRIM(ADJUSTL(binary_restart_file_name)), &
2447 136 : "Number of atomic kinds:", natomkind, &
2448 136 : "Number of atoms:", natom, &
2449 136 : "Number of cores (only core-shell model):", ncore, &
2450 136 : "Number of shells (only core-shell model):", nshell, &
2451 136 : "Number of molecule kinds:", nmoleculekind, &
2452 272 : "Number of molecules", nmolecule
2453 : END IF
2454 :
2455 136 : n_int_size = n_int_size + 6
2456 : END IF
2457 :
2458 : WRITE (UNIT=output_unit, IOSTAT=istat) &
2459 136 : natomkind, natom, ncore, nshell, nmoleculekind, nmolecule
2460 136 : IF (istat /= 0) THEN
2461 : CALL stop_write("natomkind,natom,ncore,nshell,nmoleculekind,nmolecule "// &
2462 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2463 0 : output_unit)
2464 : END IF
2465 :
2466 : ! Write atomic kind names
2467 408 : DO ikind = 1, natomkind
2468 272 : WRITE (UNIT=output_unit, IOSTAT=istat) atomic_kinds%els(ikind)%name
2469 272 : IF (istat /= 0) CALL stop_write("atomic_kinds%els(ikind)%name "// &
2470 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2471 0 : output_unit)
2472 408 : n_char_size = n_char_size + LEN(atomic_kinds%els(ikind)%name)
2473 : END DO
2474 :
2475 : ! Write atomic kind numbers of all atoms
2476 408 : ALLOCATE (ibuf(natom))
2477 13192 : DO iatom = 1, natom
2478 13192 : ibuf(iatom) = particles%els(iatom)%atomic_kind%kind_number
2479 : END DO
2480 136 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:natom)
2481 136 : IF (istat /= 0) CALL stop_write("ibuf(1:natom) -> atomic kind numbers "// &
2482 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2483 0 : output_unit)
2484 136 : n_int_size = n_int_size + natom
2485 : ! Write atomic coordinates
2486 408 : ALLOCATE (rbuf(3, natom))
2487 13192 : DO iatom = 1, natom
2488 52360 : rbuf(1:3, iatom) = particles%els(iatom)%r(1:3)
2489 : END DO
2490 136 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:natom)
2491 136 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:natom) -> atomic coordinates "// &
2492 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2493 0 : output_unit)
2494 136 : n_dp_size = n_dp_size + 3*natom
2495 136 : DEALLOCATE (rbuf)
2496 :
2497 : ! Write molecule information if available
2498 136 : IF (nmolecule > 0) THEN
2499 : ! Write molecule kind names
2500 272 : DO ikind = 1, nmoleculekind
2501 136 : WRITE (UNIT=output_unit, IOSTAT=istat) molecule_kinds%els(ikind)%name
2502 136 : IF (istat /= 0) CALL stop_write("molecule_kinds%els(ikind)%name "// &
2503 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2504 0 : output_unit)
2505 272 : n_char_size = n_char_size + LEN(molecule_kinds%els(ikind)%name)
2506 : END DO
2507 : ! Write molecule (kind) index numbers for all atoms
2508 136 : ibuf(:) = 0
2509 272 : ALLOCATE (imol(natom))
2510 136 : imol(:) = 0
2511 1224 : DO imolecule = 1, nmolecule
2512 1088 : ikind = molecules%els(imolecule)%molecule_kind%kind_number
2513 14144 : DO iatom = molecules%els(imolecule)%first_atom, &
2514 1224 : molecules%els(imolecule)%last_atom
2515 13056 : ibuf(iatom) = ikind
2516 14144 : imol(iatom) = imolecule
2517 : END DO
2518 : END DO
2519 : ! Write molecule kind index number for each atom
2520 136 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:natom)
2521 136 : IF (istat /= 0) CALL stop_write("ibuf(1:natom) -> molecule kind index numbers "// &
2522 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2523 0 : output_unit)
2524 136 : n_int_size = n_int_size + natom
2525 : ! Write molecule index number for each atom
2526 136 : WRITE (UNIT=output_unit, IOSTAT=istat) imol(1:natom)
2527 136 : IF (istat /= 0) CALL stop_write("imol(1:natom) -> molecule index numbers "// &
2528 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2529 0 : output_unit)
2530 136 : n_int_size = n_int_size + natom
2531 136 : DEALLOCATE (imol)
2532 : END IF ! molecules
2533 :
2534 136 : DEALLOCATE (ibuf)
2535 :
2536 : ! Core-shell model only
2537 136 : section_label = "SHELL COORDINATES"
2538 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nshell
2539 136 : IF (istat /= 0) CALL stop_write("section_label, nshell "// &
2540 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2541 0 : output_unit)
2542 136 : n_char_size = n_char_size + LEN(section_label)
2543 136 : n_int_size = n_int_size + 1
2544 136 : IF (nshell > 0) THEN
2545 : ! Write shell coordinates
2546 168 : ALLOCATE (rbuf(3, nshell))
2547 5432 : DO ishell = 1, nshell
2548 21560 : rbuf(1:3, ishell) = shell_particles%els(ishell)%r(1:3)
2549 : END DO
2550 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:nshell)
2551 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:nshell) -> shell coordinates "// &
2552 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2553 0 : output_unit)
2554 56 : n_dp_size = n_dp_size + 3*nshell
2555 56 : DEALLOCATE (rbuf)
2556 : ! Write atomic indices, i.e. number of the atom the shell belongs to
2557 168 : ALLOCATE (ibuf(nshell))
2558 5432 : DO ishell = 1, nshell
2559 5432 : ibuf(ishell) = shell_particles%els(ishell)%atom_index
2560 : END DO
2561 56 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:nshell)
2562 56 : IF (istat /= 0) CALL stop_write("ibuf(1:nshell) -> atomic indices "// &
2563 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2564 0 : output_unit)
2565 56 : n_int_size = n_int_size + nshell
2566 56 : DEALLOCATE (ibuf)
2567 : END IF
2568 :
2569 136 : section_label = "CORE COORDINATES"
2570 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, ncore
2571 136 : IF (istat /= 0) CALL stop_write("section_label, ncore "// &
2572 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2573 0 : output_unit)
2574 136 : n_char_size = n_char_size + LEN(section_label)
2575 136 : n_int_size = n_int_size + 1
2576 136 : IF (ncore > 0) THEN
2577 : ! Write core coordinates
2578 168 : ALLOCATE (rbuf(3, ncore))
2579 5432 : DO icore = 1, ncore
2580 21560 : rbuf(1:3, icore) = core_particles%els(icore)%r(1:3)
2581 : END DO
2582 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:ncore)
2583 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:ncore) -> core coordinates "// &
2584 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2585 0 : output_unit)
2586 56 : n_dp_size = n_dp_size + 3*ncore
2587 56 : DEALLOCATE (rbuf)
2588 : ! Write atomic indices, i.e. number of the atom the core belongs to
2589 168 : ALLOCATE (ibuf(ncore))
2590 5432 : DO icore = 1, ncore
2591 5432 : ibuf(icore) = core_particles%els(icore)%atom_index
2592 : END DO
2593 56 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:ncore)
2594 56 : IF (istat /= 0) CALL stop_write("ibuf(1:ncore) -> atomic indices "// &
2595 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2596 0 : output_unit)
2597 56 : n_int_size = n_int_size + ncore
2598 56 : DEALLOCATE (ibuf)
2599 : END IF
2600 : END IF ! ionode only
2601 :
2602 : ! Thermostat information
2603 :
2604 : ! Particle thermostats
2605 272 : section_label = "PARTICLE THERMOSTATS"
2606 272 : IF (ASSOCIATED(thermostat_part)) THEN
2607 : ! Nose-Hoover thermostats
2608 176 : IF (thermostat_part%type_of_thermostat == do_thermo_nose) THEN
2609 176 : nhc => thermostat_part%nhc
2610 : CALL write_binary_thermostats_nose(nhc, output_unit, log_unit, section_label, &
2611 : n_char_size, n_dp_size, n_int_size, &
2612 176 : print_info, para_env)
2613 : END IF
2614 : ELSE
2615 96 : nhc_size = 0
2616 96 : IF (output_unit > 0) THEN
2617 48 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nhc_size
2618 48 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", nhc_size "// &
2619 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2620 0 : output_unit)
2621 : END IF
2622 96 : n_char_size = n_char_size + LEN(section_label)
2623 96 : n_int_size = n_int_size + 1
2624 96 : IF (output_unit > 0 .AND. log_unit > 0) THEN ! only ionode
2625 48 : IF (print_info) THEN
2626 : WRITE (UNIT=log_unit, FMT="(T3,A,T71,I10)") &
2627 48 : "NHC size ("//TRIM(ADJUSTL(section_label))//")", nhc_size
2628 : END IF
2629 : END IF
2630 : END IF
2631 :
2632 : ! Shell thermostats (only for core-shell models)
2633 272 : section_label = "SHELL THERMOSTATS"
2634 272 : IF (ASSOCIATED(thermostat_shell)) THEN
2635 : ! Nose-Hoover thermostats
2636 24 : IF (thermostat_shell%type_of_thermostat == do_thermo_nose) THEN
2637 24 : nhc => thermostat_shell%nhc
2638 : CALL write_binary_thermostats_nose(nhc, output_unit, log_unit, section_label, &
2639 : n_char_size, n_dp_size, n_int_size, &
2640 24 : print_info, para_env)
2641 : END IF
2642 : ELSE
2643 248 : nhc_size = 0
2644 248 : IF (output_unit > 0) THEN
2645 124 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nhc_size
2646 124 : IF (istat /= 0) CALL stop_write("nhc_size "// &
2647 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2648 0 : output_unit)
2649 : END IF
2650 248 : n_char_size = n_char_size + LEN(section_label)
2651 248 : n_int_size = n_int_size + 1
2652 248 : IF (output_unit > 0 .AND. log_unit > 0) THEN ! only ionode
2653 124 : IF (print_info) THEN
2654 : WRITE (UNIT=log_unit, FMT="(T3,A,T71,I10)") &
2655 124 : "NHC size ("//TRIM(ADJUSTL(section_label))//")", nhc_size
2656 : END IF
2657 : END IF
2658 : END IF
2659 :
2660 : ! Particle velocities
2661 :
2662 272 : IF (output_unit > 0) THEN ! only ionode
2663 : ! Write particle velocities if needed
2664 136 : section_label = "VELOCITIES"
2665 : IF (output_unit > 0) THEN
2666 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, MERGE(natom, 0, write_velocities)
2667 136 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", write_velocities "// &
2668 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2669 0 : output_unit)
2670 : END IF
2671 136 : n_char_size = n_char_size + LEN(section_label)
2672 136 : n_int_size = n_int_size + 1
2673 136 : IF (print_info .AND. log_unit > 0) THEN
2674 : WRITE (UNIT=log_unit, FMT="(T3,A,T78,A3)") &
2675 136 : "Write "//TRIM(ADJUSTL(section_label))//" section", MERGE("YES", " NO", write_velocities)
2676 : END IF
2677 136 : IF (write_velocities) THEN
2678 408 : ALLOCATE (rbuf(3, natom))
2679 : ! Write atomic velocities
2680 13192 : DO iatom = 1, natom
2681 52360 : rbuf(1:3, iatom) = particles%els(iatom)%v(1:3)
2682 : END DO
2683 136 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:natom)
2684 136 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:natom) -> atomic velocities "// &
2685 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2686 0 : output_unit)
2687 136 : n_dp_size = n_dp_size + 3*natom
2688 136 : DEALLOCATE (rbuf)
2689 : END IF
2690 : ! Write shell velocities
2691 136 : section_label = "SHELL VELOCITIES"
2692 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, MERGE(nshell, 0, write_velocities)
2693 136 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", write_velocities "// &
2694 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2695 0 : output_unit)
2696 : n_char_size = n_char_size + LEN(section_label)
2697 : n_int_size = n_int_size + 1
2698 136 : IF (print_info .AND. log_unit > 0) THEN
2699 : WRITE (UNIT=log_unit, FMT="(T3,A,T78,A3)") &
2700 136 : "Write "//TRIM(ADJUSTL(section_label))//" section", MERGE("YES", " NO", write_velocities)
2701 : END IF
2702 136 : IF (nshell > 0) THEN
2703 56 : IF (write_velocities) THEN
2704 168 : ALLOCATE (rbuf(3, nshell))
2705 5432 : DO ishell = 1, nshell
2706 21560 : rbuf(1:3, ishell) = shell_particles%els(ishell)%v(1:3)
2707 : END DO
2708 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:nshell)
2709 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:nshell) -> shell velocities "// &
2710 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2711 0 : output_unit)
2712 56 : n_dp_size = n_dp_size + 3*nshell
2713 56 : DEALLOCATE (rbuf)
2714 : END IF
2715 : END IF
2716 : ! Write core velocities
2717 136 : section_label = "CORE VELOCITIES"
2718 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, MERGE(ncore, 0, write_velocities)
2719 136 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", write_velocities "// &
2720 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2721 0 : output_unit)
2722 136 : n_char_size = n_char_size + LEN(section_label)
2723 136 : n_int_size = n_int_size + 1
2724 136 : IF (print_info .AND. log_unit > 0) THEN
2725 : WRITE (UNIT=log_unit, FMT="(T3,A,T78,A3)") &
2726 136 : "Write "//TRIM(ADJUSTL(section_label))//" section", MERGE("YES", " NO", write_velocities)
2727 : END IF
2728 136 : IF (ncore > 0) THEN
2729 56 : IF (write_velocities) THEN
2730 168 : ALLOCATE (rbuf(3, ncore))
2731 5432 : DO icore = 1, ncore
2732 21560 : rbuf(1:3, icore) = core_particles%els(icore)%v(1:3)
2733 : END DO
2734 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:ncore)
2735 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:ncore) -> core velocities "// &
2736 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2737 0 : output_unit)
2738 56 : n_dp_size = n_dp_size + 3*ncore
2739 56 : DEALLOCATE (rbuf)
2740 : END IF
2741 : END IF
2742 : END IF ! ionode only
2743 :
2744 : ! Optionally, print a small I/O statistics
2745 272 : IF (output_unit > 0) THEN ! only ionode
2746 136 : IF (print_info .AND. log_unit > 0) THEN
2747 : WRITE (UNIT=log_unit, FMT="(/,(T2,I10,1X,I0,A,T68,I10,A))") &
2748 136 : n_char_size, int_size, "-byte characters written", n_char_size*int_size/1024, " KB", &
2749 136 : n_dp_size, dp_size, "-byte floating point numbers written", n_dp_size*dp_size/1024, " KB", &
2750 272 : n_int_size, int_size, "-byte integer numbers written", n_int_size*int_size/1024, " KB"
2751 : WRITE (UNIT=log_unit, FMT="(/,T2,A)") &
2752 136 : "Binary restart file "//TRIM(ADJUSTL(binary_restart_file_name))//" written"
2753 : END IF
2754 : END IF ! ionode only
2755 :
2756 272 : CALL timestop(handle)
2757 :
2758 : END SUBROUTINE write_binary_restart
2759 :
2760 : ! **************************************************************************************************
2761 : !> \brief Write an input section for Nose thermostats to an external file in
2762 : !> binary format
2763 : !> \param nhc ...
2764 : !> \param output_unit binary file to write to
2765 : !> \param log_unit unit for logging debug information
2766 : !> \param section_label ...
2767 : !> \param n_char_size ...
2768 : !> \param n_dp_size ...
2769 : !> \param n_int_size ...
2770 : !> \param print_info ...
2771 : !> \param para_env ...
2772 : !> \par History
2773 : !> - Creation (23.03.2011,MK)
2774 : !> \author Matthias Krack (MK)
2775 : !> \version 1.0
2776 : ! **************************************************************************************************
2777 200 : SUBROUTINE write_binary_thermostats_nose(nhc, output_unit, log_unit, section_label, &
2778 : n_char_size, n_dp_size, n_int_size, &
2779 : print_info, para_env)
2780 :
2781 : TYPE(lnhc_parameters_type), POINTER :: nhc
2782 : INTEGER, INTENT(IN) :: output_unit, log_unit
2783 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section_label
2784 : INTEGER, INTENT(INOUT) :: n_char_size, n_dp_size, n_int_size
2785 : LOGICAL, INTENT(IN) :: print_info
2786 : TYPE(mp_para_env_type), POINTER :: para_env
2787 :
2788 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_binary_thermostats_nose'
2789 :
2790 : INTEGER :: handle, istat, nhc_size
2791 200 : REAL(KIND=dp), DIMENSION(:), POINTER :: eta, fnhc, mnhc, veta
2792 :
2793 200 : CALL timeset(routineN, handle)
2794 :
2795 200 : NULLIFY (eta)
2796 200 : NULLIFY (fnhc)
2797 200 : NULLIFY (mnhc)
2798 200 : NULLIFY (veta)
2799 :
2800 200 : CALL collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
2801 :
2802 200 : nhc_size = SIZE(eta)
2803 :
2804 200 : IF (output_unit > 0) THEN ! only ionode
2805 100 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nhc_size
2806 100 : IF (istat /= 0) CALL stop_write("nhc_size "// &
2807 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2808 0 : output_unit)
2809 100 : n_char_size = n_char_size + LEN(section_label)
2810 100 : n_int_size = n_int_size + 1
2811 100 : IF (print_info .AND. log_unit > 0) THEN
2812 : WRITE (UNIT=log_unit, FMT="(T3,A,T71,I10)") &
2813 100 : "NHC size ("//TRIM(ADJUSTL(section_label))//")", nhc_size
2814 : END IF
2815 : ! eta
2816 86500 : WRITE (UNIT=output_unit, IOSTAT=istat) eta(1:nhc_size)
2817 100 : IF (istat /= 0) CALL stop_write("eta(1:nhc_size) "// &
2818 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2819 0 : output_unit)
2820 100 : n_dp_size = n_dp_size + nhc_size
2821 : END IF ! ionode only
2822 :
2823 200 : DEALLOCATE (eta)
2824 :
2825 : ! veta
2826 200 : IF (output_unit > 0) THEN ! only ionode
2827 86500 : WRITE (UNIT=output_unit, IOSTAT=istat) veta(1:nhc_size)
2828 100 : IF (istat /= 0) CALL stop_write("veta(1:nhc_size) "// &
2829 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2830 0 : output_unit)
2831 100 : n_dp_size = n_dp_size + nhc_size
2832 : END IF ! ionode only
2833 :
2834 200 : DEALLOCATE (veta)
2835 :
2836 : ! mnhc
2837 200 : IF (output_unit > 0) THEN ! only ionode
2838 86500 : WRITE (UNIT=output_unit, IOSTAT=istat) mnhc(1:nhc_size)
2839 100 : IF (istat /= 0) CALL stop_write("mnhc(1:nhc_size) "// &
2840 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2841 0 : output_unit)
2842 100 : n_dp_size = n_dp_size + nhc_size
2843 : END IF ! ionode only
2844 :
2845 200 : DEALLOCATE (mnhc)
2846 :
2847 : ! fnhc
2848 200 : IF (output_unit > 0) THEN ! only ionode
2849 86500 : WRITE (UNIT=output_unit, IOSTAT=istat) fnhc(1:nhc_size)
2850 100 : IF (istat /= 0) CALL stop_write("fnhc(1:nhc_size) "// &
2851 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2852 0 : output_unit)
2853 100 : n_dp_size = n_dp_size + nhc_size
2854 : END IF ! ionode only
2855 :
2856 200 : DEALLOCATE (fnhc)
2857 :
2858 200 : CALL timestop(handle)
2859 :
2860 200 : END SUBROUTINE write_binary_thermostats_nose
2861 :
2862 : ! **************************************************************************************************
2863 : !> \brief Print an error message and stop the program execution in case of a
2864 : !> read error.
2865 : !> \param object ...
2866 : !> \param unit_number ...
2867 : !> \par History
2868 : !> - Creation (15.02.2011,MK)
2869 : !> \author Matthias Krack (MK)
2870 : !> \note
2871 : !> object : Name of the data object for which I/O operation failed
2872 : !> unit_number: Logical unit number of the file written to
2873 : ! **************************************************************************************************
2874 0 : SUBROUTINE stop_write(object, unit_number)
2875 :
2876 : CHARACTER(LEN=*), INTENT(IN) :: object
2877 : INTEGER, INTENT(IN) :: unit_number
2878 :
2879 : CHARACTER(LEN=2*default_path_length) :: message
2880 : CHARACTER(LEN=default_path_length) :: file_name
2881 : LOGICAL :: file_exists
2882 :
2883 0 : IF (unit_number >= 0) THEN
2884 0 : INQUIRE (UNIT=unit_number, EXIST=file_exists)
2885 : ELSE
2886 0 : file_exists = .FALSE.
2887 : END IF
2888 0 : IF (file_exists) THEN
2889 0 : INQUIRE (UNIT=unit_number, NAME=file_name)
2890 : WRITE (UNIT=message, FMT="(A)") &
2891 : "An error occurred writing data object <"//TRIM(ADJUSTL(object))// &
2892 0 : "> to file <"//TRIM(ADJUSTL(file_name))//">"
2893 : ELSE
2894 : WRITE (UNIT=message, FMT="(A,I0,A)") &
2895 : "Could not write data object <"//TRIM(ADJUSTL(object))// &
2896 0 : "> to logical unit ", unit_number, ". The I/O unit does not exist."
2897 : END IF
2898 :
2899 0 : CPABORT(message)
2900 :
2901 0 : END SUBROUTINE stop_write
2902 :
2903 : END MODULE input_cp2k_restarts
|
