Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 117754 : 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 58877 : CALL timeset(routineN, handle)
130 :
131 58877 : logger => cp_get_default_logger()
132 58877 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
133 :
134 58877 : NULLIFY (global_section)
135 58877 : global_section => section_vals_get_subs_vals(root_section, "GLOBAL")
136 58877 : 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 58877 : 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 14694 : sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
144 14694 : CALL section_vals_get(sections, n_repetition=nforce_eval)
145 : CALL section_vals_val_get(motion_section, "PRINT%RESTART%SPLIT_RESTART_FILE", &
146 14694 : l_val=write_binary_restart_file)
147 :
148 14694 : 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 14694 : write_binary_restart_file=write_binary_restart_file)
172 :
173 44082 : DO ikey = 1, SIZE(keys)
174 29388 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
175 14694 : 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 15918 : do_backup=(ikey == 2))
180 15918 : IF (ires > 0) THEN
181 8312 : CALL write_restart_header(ires)
182 8312 : CALL section_vals_write(root_section, unit_nr=ires, hide_root=.TRUE.)
183 : END IF
184 15918 : CALL cp_print_key_finished_output(ires, logger, motion_section, TRIM(keys(ikey)))
185 : END IF
186 : END DO
187 :
188 14694 : IF (save_mem) THEN
189 76 : CALL update_subsys_release(md_env, force_env, root_section)
190 76 : CALL update_motion_release(motion_section)
191 : END IF
192 :
193 : END IF
194 :
195 58877 : CALL timestop(handle)
196 :
197 58877 : 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 212 : 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 212 : 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 212 : CALL timeset(routineN, handle)
228 :
229 212 : NULLIFY (core_particles, my_force_env, my_force_b, particles, &
230 212 : shell_particles, subsys, work_section)
231 :
232 212 : IF (PRESENT(md_env)) THEN
233 148 : CALL get_md_env(md_env=md_env, force_env=my_force_env)
234 64 : ELSEIF (PRESENT(force_env)) THEN
235 64 : my_force_env => force_env
236 : END IF
237 :
238 212 : IF (ASSOCIATED(my_force_env)) THEN
239 212 : NULLIFY (subsys_section)
240 212 : 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 212 : (myid /= pint_run))
245 :
246 212 : force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
247 212 : CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval)
248 :
249 424 : DO iforce_eval = 1, nforce_eval
250 : subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", &
251 212 : i_rep_section=i_force_eval(iforce_eval))
252 212 : CALL section_vals_get(subsys_section, explicit=explicit)
253 212 : IF (.NOT. explicit) CYCLE ! Nothing to update...
254 :
255 212 : my_force_b => my_force_env
256 212 : IF (iforce_eval > 1) my_force_b => my_force_env%sub_force_env(iforce_eval - 1)%force_env
257 :
258 212 : CALL force_env_get(my_force_b, subsys=subsys)
259 :
260 : CALL cp_subsys_get(subsys, particles=particles, shell_particles=shell_particles, &
261 212 : core_particles=core_particles)
262 :
263 212 : work_section => section_vals_get_subs_vals(subsys_section, "COORD")
264 212 : CALL section_vals_get(work_section, explicit=explicit)
265 212 : IF (explicit) THEN
266 212 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
267 212 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
268 : END IF
269 212 : CALL section_vals_remove_values(work_section)
270 212 : IF (explicit) THEN
271 212 : CALL section_vals_val_set(work_section, "UNIT", c_val=unit_str)
272 212 : CALL section_vals_val_set(work_section, "SCALED", l_val=scale)
273 : END IF
274 :
275 212 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
276 212 : IF (.NOT. skip_vel_section) THEN
277 148 : CALL section_vals_remove_values(work_section)
278 : END IF
279 :
280 212 : 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 848 : 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 212 : DEALLOCATE (i_force_eval)
321 :
322 : END IF
323 :
324 212 : CALL timestop(handle)
325 :
326 212 : 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 212 : 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 212 : CALL timeset(routineN, handle)
345 :
346 212 : NULLIFY (work_section)
347 :
348 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%AVERAGES%RESTART_AVERAGES")
349 212 : CALL section_vals_remove_values(work_section)
350 :
351 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%COORD")
352 212 : CALL section_vals_remove_values(work_section)
353 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%VELOCITY")
354 212 : CALL section_vals_remove_values(work_section)
355 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%MASS")
356 212 : CALL section_vals_remove_values(work_section)
357 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE%FORCE")
358 212 : CALL section_vals_remove_values(work_section)
359 :
360 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%COORD")
361 212 : CALL section_vals_remove_values(work_section)
362 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%VELOCITY")
363 212 : CALL section_vals_remove_values(work_section)
364 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%MASS")
365 212 : CALL section_vals_remove_values(work_section)
366 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%NOSE%FORCE")
367 212 : CALL section_vals_remove_values(work_section)
368 :
369 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%COORD")
370 212 : CALL section_vals_remove_values(work_section)
371 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%VELOCITY")
372 212 : CALL section_vals_remove_values(work_section)
373 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%MASS")
374 212 : CALL section_vals_remove_values(work_section)
375 212 : work_section => section_vals_get_subs_vals(motion_section, "MD%SHELL%THERMOSTAT%NOSE%FORCE")
376 212 : CALL section_vals_remove_values(work_section)
377 :
378 212 : CALL timestop(handle)
379 :
380 212 : 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 14694 : 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 14694 : CALL timeset(routineN, handle)
421 :
422 14694 : NULLIFY (logger, motion_section, my_force_env)
423 :
424 14694 : IF (PRESENT(save_mem)) THEN
425 14694 : my_save_mem = save_mem
426 : ELSE
427 0 : my_save_mem = .FALSE.
428 : END IF
429 :
430 14694 : IF (PRESENT(write_binary_restart_file)) THEN
431 14694 : my_write_binary_restart_file = write_binary_restart_file
432 : ELSE
433 0 : my_write_binary_restart_file = .FALSE.
434 : END IF
435 :
436 14694 : logger => cp_get_default_logger()
437 :
438 : ! Can handle md_env or force_env
439 14694 : lcond = PRESENT(md_env) .OR. PRESENT(force_env) .OR. PRESENT(pint_env) .OR. PRESENT(helium_env)
440 : IF (lcond) THEN
441 14568 : IF (PRESENT(md_env)) THEN
442 5508 : CALL get_md_env(md_env=md_env, force_env=my_force_env)
443 9060 : ELSE IF (PRESENT(force_env)) THEN
444 8622 : my_force_env => force_env
445 : END IF
446 : ! The real restart setting...
447 14568 : 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 29026 : write_binary_restart_file=my_write_binary_restart_file)
458 : ! Update one force_env_section per time..
459 14568 : IF (ASSOCIATED(my_force_env)) THEN
460 14130 : do_respa = .FALSE.
461 : ! Do respa only in case of RESPA MD
462 14130 : IF (PRESENT(md_env)) THEN
463 5508 : CALL get_md_env(md_env=md_env, simpar=simpar)
464 5508 : 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 14130 : respa=do_respa)
473 :
474 : END IF
475 : END IF
476 :
477 14694 : CALL timestop(handle)
478 :
479 14694 : 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 131112 : 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 14568 : INTEGER, DIMENSION(:), POINTER :: walkers_status
517 : INTEGER, POINTER :: itimes
518 : LOGICAL :: my_save_mem, my_write_binary_restart_file
519 14568 : 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 14568 : 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 14568 : CALL timeset(routineN, handle)
534 14568 : NULLIFY (logger, thermostat_part, thermostat_baro, npt, para_env, nhc, &
535 14568 : work_section, thermostat_shell, t, averages, constant, &
536 14568 : walkers_status, itimes, meta_env, simpar)
537 14568 : NULLIFY (particles)
538 14568 : NULLIFY (subsys)
539 14568 : 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 5508 : para_env=para_env)
551 : ELSE
552 9060 : IF (ASSOCIATED(force_env)) THEN
553 8622 : para_env => force_env%para_env
554 438 : ELSEIF (PRESENT(pint_env)) THEN
555 400 : para_env => pint_env%logger%para_env
556 38 : ELSEIF (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 14568 : IF (ASSOCIATED(force_env)) THEN
567 14130 : meta_env => force_env%meta_env
568 : END IF
569 :
570 : IF (PRESENT(save_mem)) THEN
571 14568 : my_save_mem = save_mem
572 : ELSE
573 : my_save_mem = .FALSE.
574 : END IF
575 :
576 14568 : IF (PRESENT(write_binary_restart_file)) THEN
577 14568 : my_write_binary_restart_file = write_binary_restart_file
578 : ELSE
579 : my_write_binary_restart_file = .FALSE.
580 : END IF
581 :
582 14568 : CALL timeset(routineN//"_COUNTERS", handle2)
583 14568 : IF (ASSOCIATED(itimes)) THEN
584 5508 : IF (itimes >= 0) THEN
585 5508 : CALL section_vals_val_set(motion_section, "MD%STEP_START_VAL", i_val=itimes)
586 5508 : CPASSERT(ASSOCIATED(t))
587 5508 : CALL section_vals_val_set(motion_section, "MD%TIME_START_VAL", r_val=t)
588 : END IF
589 : END IF
590 14568 : IF (ASSOCIATED(constant)) THEN
591 5508 : CALL section_vals_val_set(motion_section, "MD%ECONS_START_VAL", r_val=constant)
592 : END IF
593 14568 : CALL timestop(handle2)
594 : ! AVERAGES
595 14568 : CALL timeset(routineN//"_AVERAGES", handle2)
596 14568 : IF (ASSOCIATED(averages)) THEN
597 5508 : IF ((averages%do_averages) .AND. (averages%itimes_start /= -1)) THEN
598 5500 : work_section => section_vals_get_subs_vals(motion_section, "MD%AVERAGES")
599 5500 : CALL section_vals_val_set(work_section, "_SECTION_PARAMETERS_", l_val=averages%do_averages)
600 5500 : work_section => section_vals_get_subs_vals(motion_section, "MD%AVERAGES%RESTART_AVERAGES")
601 5500 : CALL section_vals_val_set(work_section, "ITIMES_START", i_val=averages%itimes_start)
602 5500 : CALL section_vals_val_set(work_section, "AVECPU", r_val=averages%avecpu)
603 5500 : CALL section_vals_val_set(work_section, "AVEHUGONIOT", r_val=averages%avehugoniot)
604 5500 : CALL section_vals_val_set(work_section, "AVETEMP_BARO", r_val=averages%avetemp_baro)
605 5500 : CALL section_vals_val_set(work_section, "AVEPOT", r_val=averages%avepot)
606 5500 : CALL section_vals_val_set(work_section, "AVEKIN", r_val=averages%avekin)
607 5500 : CALL section_vals_val_set(work_section, "AVETEMP", r_val=averages%avetemp)
608 5500 : CALL section_vals_val_set(work_section, "AVEKIN_QM", r_val=averages%avekin_qm)
609 5500 : CALL section_vals_val_set(work_section, "AVETEMP_QM", r_val=averages%avetemp_qm)
610 5500 : CALL section_vals_val_set(work_section, "AVEVOL", r_val=averages%avevol)
611 5500 : CALL section_vals_val_set(work_section, "AVECELL_A", r_val=averages%aveca)
612 5500 : CALL section_vals_val_set(work_section, "AVECELL_B", r_val=averages%avecb)
613 5500 : CALL section_vals_val_set(work_section, "AVECELL_C", r_val=averages%avecc)
614 5500 : CALL section_vals_val_set(work_section, "AVEALPHA", r_val=averages%aveal)
615 5500 : CALL section_vals_val_set(work_section, "AVEBETA", r_val=averages%avebe)
616 5500 : CALL section_vals_val_set(work_section, "AVEGAMMA", r_val=averages%avega)
617 5500 : CALL section_vals_val_set(work_section, "AVE_ECONS", r_val=averages%econs)
618 5500 : CALL section_vals_val_set(work_section, "AVE_PRESS", r_val=averages%avepress)
619 5500 : CALL section_vals_val_set(work_section, "AVE_PXX", r_val=averages%avepxx)
620 : ! Virial averages
621 5500 : IF (ASSOCIATED(averages%virial)) THEN
622 20 : ALLOCATE (buffer(9))
623 200 : buffer = RESHAPE(averages%virial%pv_total, (/9/))
624 20 : CALL section_vals_val_set(work_section, "AVE_PV_TOT", r_vals_ptr=buffer)
625 :
626 20 : ALLOCATE (buffer(9))
627 200 : buffer = RESHAPE(averages%virial%pv_virial, (/9/))
628 20 : CALL section_vals_val_set(work_section, "AVE_PV_VIR", r_vals_ptr=buffer)
629 :
630 20 : ALLOCATE (buffer(9))
631 200 : buffer = RESHAPE(averages%virial%pv_kinetic, (/9/))
632 20 : CALL section_vals_val_set(work_section, "AVE_PV_KIN", r_vals_ptr=buffer)
633 :
634 20 : ALLOCATE (buffer(9))
635 200 : buffer = RESHAPE(averages%virial%pv_constraint, (/9/))
636 20 : CALL section_vals_val_set(work_section, "AVE_PV_CNSTR", r_vals_ptr=buffer)
637 :
638 20 : ALLOCATE (buffer(9))
639 200 : buffer = RESHAPE(averages%virial%pv_xc, (/9/))
640 20 : CALL section_vals_val_set(work_section, "AVE_PV_XC", r_vals_ptr=buffer)
641 :
642 20 : ALLOCATE (buffer(9))
643 200 : buffer = RESHAPE(averages%virial%pv_fock_4c, (/9/))
644 20 : CALL section_vals_val_set(work_section, "AVE_PV_FOCK_4C", r_vals_ptr=buffer)
645 : END IF
646 : ! Colvars averages
647 5500 : 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 5500 : 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 14568 : CALL timestop(handle2)
660 :
661 : ! SAVE THERMOSTAT target TEMPERATURE when doing TEMPERATURE_ANNEALING
662 14568 : IF (PRESENT(md_env)) THEN
663 5508 : IF (ASSOCIATED(simpar)) THEN
664 5508 : 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 14568 : CALL timeset(routineN//"_THERMOSTAT_PARTICLES", handle2)
672 14568 : IF (ASSOCIATED(thermostat_part)) THEN
673 1058 : IF (thermostat_part%type_of_thermostat == do_thermo_nose) THEN
674 : ! Restart of Nose-Hoover Thermostat for Particles
675 746 : IF (.NOT. my_write_binary_restart_file) THEN
676 658 : nhc => thermostat_part%nhc
677 658 : CALL collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
678 658 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%NOSE")
679 658 : CALL set_template_restart(work_section, eta, veta, fnhc, mnhc)
680 : END IF
681 312 : ELSE IF (thermostat_part%type_of_thermostat == do_thermo_csvr) THEN
682 : ! Restart of CSVR Thermostat for Particles
683 292 : work_section => section_vals_get_subs_vals(motion_section, "MD%THERMOSTAT%CSVR")
684 292 : 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 14568 : CALL timestop(handle2)
696 :
697 : ! BAROSTAT - THERMOSTAT
698 14568 : CALL timeset(routineN//"_BAROSTAT", handle2)
699 14568 : IF (ASSOCIATED(thermostat_baro)) THEN
700 334 : 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 756 : ALLOCATE (veta(tot_nhcneed*nhc_len))
707 756 : ALLOCATE (fnhc(tot_nhcneed*nhc_len))
708 756 : 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 82 : ELSE IF (thermostat_baro%type_of_thermostat == do_thermo_csvr) THEN
722 : ! Restart of CSVR Thermostat for Barostat
723 82 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT%THERMOSTAT%CSVR")
724 82 : CALL dump_csvr_restart_info(thermostat_baro%csvr, para_env, work_section)
725 : END IF
726 : END IF
727 14568 : CALL timestop(handle2)
728 :
729 : ! BAROSTAT
730 14568 : CALL timeset(routineN//"_NPT", handle2)
731 14568 : IF (ASSOCIATED(npt)) THEN
732 1188 : ALLOCATE (veta(SIZE(npt, 1)*SIZE(npt, 2)))
733 1188 : ALLOCATE (mnhc(SIZE(npt, 1)*SIZE(npt, 2)))
734 396 : counter = 0
735 1116 : DO i = 1, SIZE(npt, 1)
736 2808 : DO j = 1, SIZE(npt, 2)
737 1692 : counter = counter + 1
738 1692 : veta(counter) = npt(i, j)%v
739 2412 : mnhc(counter) = npt(i, j)%mass
740 : END DO
741 : END DO
742 396 : work_section => section_vals_get_subs_vals(motion_section, "MD%BAROSTAT")
743 396 : CALL set_template_restart(work_section, veta=veta, mnhc=mnhc)
744 : END IF
745 14568 : CALL timestop(handle2)
746 :
747 : ! SHELL THERMOSTAT
748 14568 : CALL timeset(routineN//"_THERMOSTAT_SHELL", handle2)
749 14568 : 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 14568 : CALL timestop(handle2)
765 :
766 14568 : CALL timeset(routineN//"_META", handle2)
767 14568 : IF (ASSOCIATED(meta_env)) THEN
768 : CALL section_vals_val_set(meta_env%metadyn_section, "STEP_START_VAL", &
769 982 : i_val=meta_env%n_steps)
770 : CALL section_vals_val_set(meta_env%metadyn_section, "NHILLS_START_VAL", &
771 982 : 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 982 : r_val=meta_env%hills_env%min_disp)
775 : CALL section_vals_val_set(meta_env%metadyn_section, "OLD_HILL_NUMBER", &
776 982 : i_val=meta_env%hills_env%old_hill_number)
777 : CALL section_vals_val_set(meta_env%metadyn_section, "OLD_HILL_STEP", &
778 982 : i_val=meta_env%hills_env%old_hill_step)
779 : !RG Adaptive hills
780 982 : IF (meta_env%do_hills .AND. meta_env%hills_env%n_hills /= 0) THEN
781 782 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_POS")
782 782 : CALL meta_hills_val_set_ss(work_section, meta_env)
783 782 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_SCALE")
784 782 : CALL meta_hills_val_set_ds(work_section, meta_env)
785 782 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "SPAWNED_HILLS_HEIGHT")
786 782 : CALL meta_hills_val_set_ww(work_section, meta_env)
787 782 : 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 982 : IF (meta_env%extended_lagrange) THEN
793 : CALL section_vals_val_set(meta_env%metadyn_section, "COLVAR_AVG_TEMPERATURE_RESTART", &
794 132 : r_val=meta_env%avg_temp)
795 132 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
796 292 : 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 292 : i_rep_val=irep)
799 : END DO
800 132 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
801 292 : 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 292 : i_rep_val=irep)
804 : END DO
805 :
806 132 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS")
807 292 : 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 292 : i_rep_val=irep)
810 : END DO
811 132 : work_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_FS")
812 292 : 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 292 : i_rep_val=irep)
815 : END DO
816 :
817 : END IF
818 : ! Multiple Walkers
819 982 : 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 14568 : CALL timestop(handle2)
827 14568 : CALL timeset(routineN//"_NEB", handle2)
828 14568 : 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 14568 : CALL timestop(handle2)
876 :
877 14568 : IF (PRESENT(pint_env)) THEN
878 : ! Update PINT section
879 400 : CALL update_motion_pint(motion_section, pint_env)
880 : END IF
881 :
882 14568 : IF (PRESENT(helium_env)) THEN
883 : ! Update HELIUM section
884 110 : CALL update_motion_helium(helium_env)
885 : END IF
886 :
887 14568 : CALL timestop(handle)
888 :
889 14568 : 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 400 : 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 400 : REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals
910 : TYPE(section_vals_type), POINTER :: pint_section, tmpsec
911 :
912 400 : CALL timeset(routineN, handle)
913 :
914 400 : pint_section => section_vals_get_subs_vals(motion_section, "PINT")
915 400 : 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 400 : NULLIFY (tmpsec)
921 400 : tmpsec => section_vals_get_subs_vals(pint_section, "BEADS")
922 400 : CALL section_vals_get(tmpsec, explicit=explicit)
923 400 : IF (.NOT. explicit) THEN
924 44 : CALL section_vals_add_values(tmpsec)
925 : END IF
926 :
927 : ! update bead coordinates in the global input structure
928 400 : NULLIFY (r_vals)
929 1200 : ALLOCATE (r_vals(pint_env%p*pint_env%ndim))
930 :
931 400 : i = 1
932 400 : CALL pint_u2x(pint_env)
933 96160 : DO iatom = 1, pint_env%ndim
934 491872 : DO ibead = 1, pint_env%p
935 395712 : r_vals(i) = pint_env%x(ibead, iatom)
936 491472 : i = i + 1
937 : END DO
938 : END DO
939 : CALL section_vals_val_set(pint_section, "BEADS%COORD%_DEFAULT_KEYWORD_", &
940 400 : r_vals_ptr=r_vals)
941 :
942 : ! update bead velocities in the global input structure
943 400 : NULLIFY (r_vals)
944 1200 : ALLOCATE (r_vals(pint_env%p*pint_env%ndim))
945 400 : i = 1
946 400 : CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
947 96160 : DO iatom = 1, pint_env%ndim
948 491872 : DO ibead = 1, pint_env%p
949 395712 : r_vals(i) = pint_env%v(ibead, iatom)
950 491472 : i = i + 1
951 : END DO
952 : END DO
953 : CALL section_vals_val_set(pint_section, "BEADS%VELOCITY%_DEFAULT_KEYWORD_", &
954 400 : r_vals_ptr=r_vals)
955 :
956 400 : 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 174 : ELSEIF (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 174 : ELSE IF (pint_env%pimd_thermostat == thermostat_pile) THEN
1005 : tmpsec => section_vals_get_subs_vals(pint_section, &
1006 102 : "PILE%RNG_INIT")
1007 102 : CALL pint_env%pile_therm%gaussian_rng_stream%dump(rng_record)
1008 102 : CALL string_to_ascii(rng_record, ascii(:, 1))
1009 : CALL section_rng_val_set(rng_section=tmpsec, nsize=1, &
1010 102 : ascii=ascii)
1011 102 : tmpsec => section_vals_get_subs_vals(pint_section, "PILE")
1012 : CALL section_vals_val_set(tmpsec, "THERMOSTAT_ENERGY", &
1013 102 : r_val=pint_env%e_pile)
1014 72 : 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 52 : 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 400 : CALL timestop(handle)
1051 :
1052 800 : 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 105 : 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 680553 : 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 .LT. 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 : 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 : 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 : 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 11040 : 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 67 : 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) &
1330 0 : CPABORT(err_str)
1331 :
1332 : ! allocate the buffer to be saved and fill it with forces
1333 : ! forces should be the same on all processors, but we don't check that here
1334 67 : NULLIFY (real_msg_gather)
1335 201 : ALLOCATE (real_msg_gather(msglen))
1336 4531 : real_msg_gather(:) = PACK(helium_env(1)%helium%force_avrg, .TRUE.)
1337 :
1338 : ! update forces in the global input structure
1339 : CALL section_vals_val_set(helium_env(1)%helium%input, &
1340 : "MOTION%PINT%HELIUM%FORCE%_DEFAULT_KEYWORD_", &
1341 67 : r_vals_ptr=real_msg_gather)
1342 :
1343 : ! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
1344 : ! assigned in section_vals_val_set - this memeory will be used later on!
1345 : ! "The val becomes the owner of the array" - from section_vals_val_set docu
1346 67 : NULLIFY (real_msg_gather)
1347 : END IF
1348 :
1349 : !
1350 : ! save the RDFs
1351 : !
1352 105 : IF (helium_env(1)%helium%rdf_present) THEN
1353 :
1354 : ! work on the temporary array so that accumulated data remains intact
1355 5010 : helium_env(1)%helium%rdf_inst(:, :) = 0.0_dp
1356 20 : DO i = 1, SIZE(helium_env)
1357 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :) + &
1358 5020 : helium_env(i)%helium%rdf_accu(:, :)
1359 : END DO
1360 :
1361 : ! average over processors / helium environments
1362 10010 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rdf_inst)
1363 10 : itmp = helium_env(1)%helium%num_env
1364 10 : invproc = 1.0_dp/REAL(itmp, dp)
1365 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*invproc
1366 :
1367 10 : nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1368 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps, dp)
1369 10 : iweight = helium_env(1)%helium%rdf_iweight
1370 : ! average over the old and the current density (observe the weights!)
1371 : helium_env(1)%helium%rdf_inst(:, :) = nsteps*helium_env(1)%helium%rdf_inst(:, :) + &
1372 5010 : iweight*helium_env(1)%helium%rdf_rstr(:, :)
1373 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps + iweight, dp)
1374 : ! update in the global input structure
1375 10 : NULLIFY (real_msg)
1376 10 : msglen = SIZE(helium_env(1)%helium%rdf_inst)
1377 30 : ALLOCATE (real_msg(msglen))
1378 2510 : real_msg(:) = PACK(helium_env(1)%helium%rdf_inst, .TRUE.)
1379 : CALL section_vals_val_set( &
1380 : helium_env(1)%helium%input, &
1381 : "MOTION%PINT%HELIUM%AVERAGES%RDF", &
1382 10 : r_vals_ptr=real_msg)
1383 10 : NULLIFY (real_msg)
1384 :
1385 : END IF
1386 :
1387 : !
1388 : ! save the densities
1389 : !
1390 105 : IF (helium_env(1)%helium%rho_present) THEN
1391 :
1392 : ! work on the temporary array so that accumulated data remains intact
1393 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = 0.0_dp
1394 20 : DO i = 1, SIZE(helium_env)
1395 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :) + &
1396 21120 : helium_env(i)%helium%rho_accu(:, :, :, :)
1397 : END DO
1398 :
1399 : ! average over processors / helium environments
1400 42210 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rho_inst)
1401 10 : itmp = helium_env(1)%helium%num_env
1402 10 : invproc = 1.0_dp/REAL(itmp, dp)
1403 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*invproc
1404 :
1405 10 : nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1406 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/REAL(nsteps, dp)
1407 10 : iweight = helium_env(1)%helium%averages_iweight
1408 : ! average over the old and the current density (observe the weights!)
1409 : helium_env(1)%helium%rho_inst(:, :, :, :) = nsteps*helium_env(1)%helium%rho_inst(:, :, :, :) + &
1410 21110 : iweight*helium_env(1)%helium%rho_rstr(:, :, :, :)
1411 21110 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/REAL(nsteps + iweight, dp)
1412 :
1413 : ! update the densities in the global input structure
1414 10 : NULLIFY (real_msg)
1415 10 : msglen = SIZE(helium_env(1)%helium%rho_inst)
1416 30 : ALLOCATE (real_msg(msglen))
1417 10010 : real_msg(:) = PACK(helium_env(1)%helium%rho_inst, .TRUE.)
1418 : CALL section_vals_val_set( &
1419 : helium_env(1)%helium%input, &
1420 : "MOTION%PINT%HELIUM%AVERAGES%RHO", &
1421 10 : r_vals_ptr=real_msg)
1422 10 : NULLIFY (real_msg)
1423 :
1424 : END IF
1425 :
1426 : END IF ! ASSOCIATED(helium_env)
1427 :
1428 110 : CALL timestop(handle)
1429 :
1430 110 : END SUBROUTINE update_motion_helium
1431 :
1432 : ! **************************************************************************************************
1433 : !> \brief routine to dump thermostat CSVR energies
1434 : !> \param thermostat_energy ...
1435 : !> \param nsize ...
1436 : !> \param work_section ...
1437 : !> \par History
1438 : !> 10.2007 created [teo]
1439 : !> \author Teodoro Laino - University of Zurich
1440 : ! **************************************************************************************************
1441 418 : SUBROUTINE dump_csvr_energy_info(thermostat_energy, nsize, work_section)
1442 :
1443 : REAL(KIND=dp), DIMENSION(:), POINTER :: thermostat_energy
1444 : INTEGER, INTENT(IN) :: nsize
1445 : TYPE(section_vals_type), POINTER :: work_section
1446 :
1447 : INTEGER :: ik, irk, Nlist
1448 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
1449 : TYPE(section_type), POINTER :: section
1450 : TYPE(val_type), POINTER :: my_val, old_val
1451 :
1452 418 : CPASSERT(ASSOCIATED(work_section))
1453 418 : CPASSERT(work_section%ref_count > 0)
1454 :
1455 418 : NULLIFY (my_val, old_val, section, vals)
1456 :
1457 418 : section => work_section%section
1458 :
1459 418 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
1460 :
1461 418 : IF (ik == -2) &
1462 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
1463 0 : "_DEFAULT_KEYWORD_")
1464 :
1465 132 : DO
1466 550 : IF (SIZE(work_section%values, 2) == 1) EXIT
1467 132 : CALL section_vals_add_values(work_section)
1468 : END DO
1469 :
1470 418 : vals => work_section%values(ik, 1)%list
1471 418 : Nlist = 0
1472 :
1473 418 : IF (ASSOCIATED(vals)) THEN
1474 286 : Nlist = cp_sll_val_get_length(vals)
1475 : END IF
1476 :
1477 22790 : DO irk = 1, nsize
1478 22372 : CALL val_create(val=my_val, r_val=thermostat_energy(irk))
1479 22372 : IF (Nlist /= 0) THEN
1480 19024 : IF (irk == 1) THEN
1481 286 : new_pos => vals
1482 : ELSE
1483 18738 : new_pos => new_pos%rest
1484 : END IF
1485 19024 : old_val => new_pos%first_el
1486 19024 : CALL val_release(old_val)
1487 19024 : new_pos%first_el => my_val
1488 : ELSE
1489 3348 : IF (irk == 1) THEN
1490 132 : NULLIFY (new_pos)
1491 132 : CALL cp_sll_val_create(new_pos, first_el=my_val)
1492 132 : vals => new_pos
1493 : ELSE
1494 3216 : NULLIFY (new_pos%rest)
1495 3216 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
1496 3216 : new_pos => new_pos%rest
1497 : END IF
1498 : END IF
1499 22790 : NULLIFY (my_val)
1500 : END DO
1501 418 : work_section%values(ik, 1)%list => vals
1502 :
1503 418 : END SUBROUTINE dump_csvr_energy_info
1504 :
1505 : ! **************************************************************************************************
1506 : !> \brief Collect all information needed to dump the restart for CSVR
1507 : !> thermostat
1508 : !> \param csvr ...
1509 : !> \param para_env ...
1510 : !> \param csvr_section ...
1511 : !> \par History
1512 : !> 10.2007 created [tlaino] - University of Zurich
1513 : !> \author Teodoro Laino
1514 : ! **************************************************************************************************
1515 398 : SUBROUTINE dump_csvr_restart_info(csvr, para_env, csvr_section)
1516 :
1517 : TYPE(csvr_system_type), POINTER :: csvr
1518 : TYPE(mp_para_env_type), POINTER :: para_env
1519 : TYPE(section_vals_type), POINTER :: csvr_section
1520 :
1521 : CHARACTER(LEN=rng_record_length) :: rng_record
1522 : INTEGER :: i, my_index
1523 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dwork
1524 : REAL(KIND=dp) :: dum
1525 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: thermo_energy
1526 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
1527 : TYPE(section_vals_type), POINTER :: work_section
1528 :
1529 : ! Thermostat Energies
1530 :
1531 1194 : ALLOCATE (work(csvr%glob_num_csvr))
1532 :
1533 1194 : ALLOCATE (thermo_energy(csvr%loc_num_csvr))
1534 8496 : DO i = 1, csvr%loc_num_csvr
1535 8496 : thermo_energy(i) = csvr%nvt(i)%thermostat_energy
1536 : END DO
1537 : CALL get_kin_energies(csvr%map_info, csvr%loc_num_csvr, &
1538 : csvr%glob_num_csvr, thermo_energy, &
1539 398 : dum, para_env, array_kin=work)
1540 398 : DEALLOCATE (thermo_energy)
1541 :
1542 : ! If check passes then let's dump the info on the restart file
1543 398 : work_section => section_vals_get_subs_vals(csvr_section, "THERMOSTAT_ENERGY")
1544 398 : CALL dump_csvr_energy_info(work, csvr%glob_num_csvr, work_section)
1545 398 : DEALLOCATE (work)
1546 :
1547 : ! Thermostat Random Number info for restart
1548 398 : work_section => section_vals_get_subs_vals(csvr_section, "RNG_INIT")
1549 1194 : ALLOCATE (dwork(rng_record_length, csvr%glob_num_csvr))
1550 6897526 : dwork = 0
1551 8496 : DO i = 1, csvr%loc_num_csvr
1552 8098 : my_index = csvr%map_info%index(i)
1553 8098 : CALL csvr%nvt(i)%gaussian_rng_stream%dump(rng_record)
1554 8496 : CALL string_to_ascii(rng_record, dwork(:, my_index))
1555 : END DO
1556 :
1557 : ! Collect data if there was no communication in this thermostat
1558 398 : IF (csvr%map_info%dis_type == do_thermo_no_communication) THEN
1559 : ! Collect data if there was no communication in this thermostat
1560 148 : CALL para_env%sum(dwork)
1561 : ELSE
1562 : ! Perform some check and collect data in case of communicating thermostats
1563 250 : CALL communication_thermo_low2(dwork, rng_record_length, csvr%glob_num_csvr, para_env)
1564 : END IF
1565 398 : CALL section_rng_val_set(rng_section=work_section, nsize=csvr%glob_num_csvr, ascii=dwork)
1566 398 : DEALLOCATE (dwork)
1567 :
1568 796 : END SUBROUTINE dump_csvr_restart_info
1569 :
1570 : ! **************************************************************************************************
1571 : !> \brief Collect all information needed to dump the restart for AD_LANGEVIN
1572 : !> thermostat
1573 : !> \param al ...
1574 : !> \param para_env ...
1575 : !> \param al_section ...
1576 : !> \par History
1577 : !> 10.2007 created [tlaino] - University of Zurich
1578 : !> \author Teodoro Laino
1579 : ! **************************************************************************************************
1580 0 : SUBROUTINE dump_al_restart_info(al, para_env, al_section)
1581 :
1582 : TYPE(al_system_type), POINTER :: al
1583 : TYPE(mp_para_env_type), POINTER :: para_env
1584 : TYPE(section_vals_type), POINTER :: al_section
1585 :
1586 : INTEGER :: i
1587 : REAL(KIND=dp) :: dum
1588 : REAL(KIND=dp), DIMENSION(:), POINTER :: t_array, work
1589 : TYPE(section_vals_type), POINTER :: work_section
1590 :
1591 : ! chi and mass
1592 :
1593 0 : ALLOCATE (work(al%glob_num_al))
1594 0 : ALLOCATE (t_array(al%loc_num_al))
1595 :
1596 : ! copy chi into temporary
1597 0 : DO i = 1, al%loc_num_al
1598 0 : t_array(i) = al%nvt(i)%chi
1599 : END DO
1600 : ! consolidate into work
1601 : CALL get_kin_energies(al%map_info, al%loc_num_al, &
1602 : al%glob_num_al, t_array, &
1603 0 : dum, para_env, array_kin=work)
1604 :
1605 : ! If check passes then let's dump the info on the restart file
1606 0 : work_section => section_vals_get_subs_vals(al_section, "CHI")
1607 0 : CALL dump_csvr_energy_info(work, al%glob_num_al, work_section)
1608 :
1609 : ! copy mass into temporary
1610 0 : DO i = 1, al%loc_num_al
1611 0 : t_array(i) = al%nvt(i)%mass
1612 : END DO
1613 : ! consolidate into work
1614 : CALL get_kin_energies(al%map_info, al%loc_num_al, &
1615 : al%glob_num_al, t_array, &
1616 0 : dum, para_env, array_kin=work)
1617 :
1618 : ! If check passes then let's dump the info on the restart file
1619 0 : work_section => section_vals_get_subs_vals(al_section, "MASS")
1620 0 : CALL dump_csvr_energy_info(work, al%glob_num_al, work_section)
1621 :
1622 0 : DEALLOCATE (t_array)
1623 0 : DEALLOCATE (work)
1624 :
1625 0 : END SUBROUTINE dump_al_restart_info
1626 :
1627 : ! **************************************************************************************************
1628 : !> \brief Collect all information needed to dump the restart for GLE
1629 : !> thermostat
1630 : !> \param gle ...
1631 : !> \param para_env ...
1632 : !> \param gle_section ...
1633 : !> \author MI
1634 : ! **************************************************************************************************
1635 20 : SUBROUTINE dump_gle_restart_info(gle, para_env, gle_section)
1636 :
1637 : TYPE(gle_type), POINTER :: gle
1638 : TYPE(mp_para_env_type), POINTER :: para_env
1639 : TYPE(section_vals_type), POINTER :: gle_section
1640 :
1641 : CHARACTER(LEN=rng_record_length) :: rng_record
1642 : INTEGER :: counter, glob_num, i, iproc, j, loc_num
1643 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dwork
1644 20 : INTEGER, DIMENSION(:), POINTER :: gle_per_proc, index
1645 : REAL(dp) :: dum
1646 20 : REAL(dp), DIMENSION(:), POINTER :: s_tmp
1647 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: thermo_energy
1648 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
1649 : TYPE(section_vals_type), POINTER :: work_section
1650 :
1651 : ! Thermostat Energies
1652 :
1653 60 : ALLOCATE (work(gle%glob_num_gle))
1654 60 : ALLOCATE (thermo_energy(gle%loc_num_gle))
1655 3260 : DO i = 1, gle%loc_num_gle
1656 3260 : thermo_energy(i) = gle%nvt(i)%thermostat_energy
1657 : END DO
1658 : CALL get_kin_energies(gle%map_info, gle%loc_num_gle, &
1659 : gle%glob_num_gle, thermo_energy, &
1660 20 : dum, para_env, array_kin=work)
1661 20 : DEALLOCATE (thermo_energy)
1662 :
1663 : ! If check passes then let's dump the info on the restart file
1664 20 : work_section => section_vals_get_subs_vals(gle_section, "THERMOSTAT_ENERGY")
1665 20 : CALL dump_csvr_energy_info(work, gle%glob_num_gle, work_section)
1666 20 : DEALLOCATE (work)
1667 :
1668 : ! Thermostat Random Number info for restart
1669 20 : work_section => section_vals_get_subs_vals(gle_section, "RNG_INIT")
1670 20 : glob_num = gle%glob_num_gle
1671 20 : loc_num = gle%loc_num_gle
1672 60 : ALLOCATE (dwork(rng_record_length, glob_num))
1673 2812340 : dwork = 0
1674 3260 : DO i = 1, loc_num
1675 3240 : j = gle%map_info%index(i)
1676 3240 : CALL gle%nvt(i)%gaussian_rng_stream%dump(rng_record)
1677 3260 : CALL string_to_ascii(rng_record, dwork(:, j))
1678 : END DO
1679 :
1680 : ! Collect data if there was no communication in this thermostat
1681 20 : IF (gle%map_info%dis_type == do_thermo_no_communication) THEN
1682 : ! Collect data if there was no communication in this thermostat
1683 20 : CALL para_env%sum(dwork)
1684 : ELSE
1685 : ! Perform some check and collect data in case of communicating thermostats
1686 0 : CALL communication_thermo_low2(dwork, rng_record_length, glob_num, para_env)
1687 : END IF
1688 20 : CALL section_rng_val_set(rng_section=work_section, nsize=glob_num, ascii=dwork)
1689 20 : DEALLOCATE (dwork)
1690 :
1691 60 : ALLOCATE (gle_per_proc(para_env%num_pe))
1692 60 : gle_per_proc(:) = 0
1693 60 : CALL para_env%allgather(gle%loc_num_gle, gle_per_proc)
1694 :
1695 : ! Thermostat S variable info for restart
1696 20 : NULLIFY (s_tmp)
1697 60 : ALLOCATE (s_tmp((gle%ndim)*gle%glob_num_gle))
1698 32420 : s_tmp = 0.0_dp
1699 :
1700 20 : NULLIFY (work, index)
1701 60 : DO iproc = 1, para_env%num_pe
1702 40 : CALL reallocate(work, 1, gle_per_proc(iproc)*(gle%ndim))
1703 40 : CALL reallocate(index, 1, gle_per_proc(iproc))
1704 40 : IF (para_env%mepos == (iproc - 1)) THEN
1705 3260 : INDEX(:) = 0
1706 20 : counter = 0
1707 120 : DO i = 1, gle%ndim
1708 16320 : DO j = 1, gle%loc_num_gle
1709 16200 : counter = counter + 1
1710 16200 : work(counter) = gle%nvt(j)%s(i)
1711 16300 : INDEX(j) = gle%map_info%index(j)
1712 : END DO
1713 : END DO
1714 : ELSE
1715 16220 : work(:) = 0.0_dp
1716 : END IF
1717 64840 : CALL para_env%bcast(work, iproc - 1)
1718 13000 : CALL para_env%bcast(index, iproc - 1)
1719 40 : counter = 0
1720 260 : DO i = 1, gle%ndim
1721 32640 : DO j = 1, gle_per_proc(iproc)
1722 32400 : counter = counter + 1
1723 32600 : s_tmp((INDEX(j) - 1)*(gle%ndim) + i) = work(counter)
1724 : END DO
1725 : END DO
1726 : END DO
1727 :
1728 20 : IF (SIZE(s_tmp) > 0) THEN
1729 20 : work_section => section_vals_get_subs_vals(gle_section, "S")
1730 20 : CALL section_vals_val_set(work_section, "_DEFAULT_KEYWORD_", r_vals_ptr=s_tmp)
1731 : ELSE
1732 0 : DEALLOCATE (s_tmp)
1733 : END IF
1734 :
1735 20 : DEALLOCATE (gle_per_proc)
1736 20 : DEALLOCATE (work)
1737 20 : DEALLOCATE (index)
1738 :
1739 40 : END SUBROUTINE dump_gle_restart_info
1740 :
1741 : ! **************************************************************************************************
1742 : !> \brief Collect all information needed to dump the restart for Nose-Hoover
1743 : !> thermostat
1744 : !> \param nhc ...
1745 : !> \param para_env ...
1746 : !> \param eta ...
1747 : !> \param veta ...
1748 : !> \param fnhc ...
1749 : !> \param mnhc ...
1750 : !> \par History
1751 : !> 10.2007 created [tlaino] - University of Zurich
1752 : !> \author Teodoro Laino
1753 : ! **************************************************************************************************
1754 982 : SUBROUTINE collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
1755 :
1756 : TYPE(lnhc_parameters_type), POINTER :: nhc
1757 : TYPE(mp_para_env_type), POINTER :: para_env
1758 : REAL(KIND=dp), DIMENSION(:), POINTER :: eta, veta, fnhc, mnhc
1759 :
1760 : INTEGER :: counter, i, iproc, j, nhc_len, num_nhc, &
1761 : numneed, tot_nhcneed
1762 982 : INTEGER, DIMENSION(:), POINTER :: index, nhc_per_proc
1763 982 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
1764 : TYPE(map_info_type), POINTER :: map_info
1765 :
1766 982 : nhc_len = SIZE(nhc%nvt, 1)
1767 982 : num_nhc = nhc%loc_num_nhc
1768 982 : numneed = num_nhc
1769 982 : map_info => nhc%map_info
1770 2946 : ALLOCATE (nhc_per_proc(para_env%num_pe))
1771 2946 : nhc_per_proc(:) = 0
1772 :
1773 2946 : CALL para_env%allgather(numneed, nhc_per_proc)
1774 982 : tot_nhcneed = nhc%glob_num_nhc
1775 :
1776 982 : NULLIFY (work, index)
1777 : !-----------------------------------------------------------------------------
1778 : !-----------------------------------------------------------------------------
1779 : ! nhc%eta
1780 : !-----------------------------------------------------------------------------
1781 2946 : ALLOCATE (eta(tot_nhcneed*nhc_len))
1782 2946 : DO iproc = 1, para_env%num_pe
1783 1964 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1784 1964 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1785 1964 : IF (para_env%mepos == (iproc - 1)) THEN
1786 68979 : INDEX(:) = 0
1787 : counter = 0
1788 4742 : DO i = 1, nhc_len
1789 257321 : DO j = 1, num_nhc
1790 252579 : counter = counter + 1
1791 252579 : work(counter) = nhc%nvt(i, j)%eta
1792 256339 : INDEX(j) = map_info%index(j)
1793 : END DO
1794 : END DO
1795 : ELSE
1796 253561 : work(:) = 0.0_dp
1797 : END IF
1798 1012280 : CALL para_env%bcast(work, iproc - 1)
1799 273952 : CALL para_env%bcast(index, iproc - 1)
1800 1964 : counter = 0
1801 10466 : DO i = 1, nhc_len
1802 514642 : DO j = 1, nhc_per_proc(iproc)
1803 505158 : counter = counter + 1
1804 512678 : eta((INDEX(j) - 1)*nhc_len + i) = work(counter)
1805 : END DO
1806 : END DO
1807 : END DO
1808 : !-----------------------------------------------------------------------------
1809 : !-----------------------------------------------------------------------------
1810 : ! nhc%veta
1811 : !-----------------------------------------------------------------------------
1812 2946 : ALLOCATE (veta(tot_nhcneed*nhc_len))
1813 2946 : DO iproc = 1, para_env%num_pe
1814 1964 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1815 1964 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1816 1964 : IF (para_env%mepos == (iproc - 1)) THEN
1817 68979 : INDEX(:) = 0
1818 : counter = 0
1819 4742 : DO i = 1, nhc_len
1820 257321 : DO j = 1, num_nhc
1821 252579 : counter = counter + 1
1822 252579 : work(counter) = nhc%nvt(i, j)%v
1823 256339 : INDEX(j) = map_info%index(j)
1824 : END DO
1825 : END DO
1826 : ELSE
1827 253561 : work(:) = 0.0_dp
1828 : END IF
1829 1012280 : CALL para_env%bcast(work, iproc - 1)
1830 273952 : CALL para_env%bcast(index, iproc - 1)
1831 1964 : counter = 0
1832 10466 : DO i = 1, nhc_len
1833 514642 : DO j = 1, nhc_per_proc(iproc)
1834 505158 : counter = counter + 1
1835 512678 : veta((INDEX(j) - 1)*nhc_len + i) = work(counter)
1836 : END DO
1837 : END DO
1838 : END DO
1839 : !-----------------------------------------------------------------------------
1840 : !-----------------------------------------------------------------------------
1841 : ! nhc%force
1842 : !-----------------------------------------------------------------------------
1843 2946 : ALLOCATE (fnhc(tot_nhcneed*nhc_len))
1844 2946 : DO iproc = 1, para_env%num_pe
1845 1964 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1846 1964 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1847 1964 : IF (para_env%mepos == (iproc - 1)) THEN
1848 68979 : INDEX(:) = 0
1849 : counter = 0
1850 4742 : DO i = 1, nhc_len
1851 257321 : DO j = 1, num_nhc
1852 252579 : counter = counter + 1
1853 252579 : work(counter) = nhc%nvt(i, j)%f
1854 256339 : INDEX(j) = map_info%index(j)
1855 : END DO
1856 : END DO
1857 : ELSE
1858 253561 : work(:) = 0.0_dp
1859 : END IF
1860 1012280 : CALL para_env%bcast(work, iproc - 1)
1861 273952 : CALL para_env%bcast(index, iproc - 1)
1862 1964 : counter = 0
1863 10466 : DO i = 1, nhc_len
1864 514642 : DO j = 1, nhc_per_proc(iproc)
1865 505158 : counter = counter + 1
1866 512678 : fnhc((INDEX(j) - 1)*nhc_len + i) = work(counter)
1867 : END DO
1868 : END DO
1869 : END DO
1870 : !-----------------------------------------------------------------------------
1871 : !-----------------------------------------------------------------------------
1872 : ! nhc%mass
1873 : !-----------------------------------------------------------------------------
1874 2946 : ALLOCATE (mnhc(tot_nhcneed*nhc_len))
1875 2946 : DO iproc = 1, para_env%num_pe
1876 1964 : CALL reallocate(work, 1, nhc_per_proc(iproc)*nhc_len)
1877 1964 : CALL reallocate(index, 1, nhc_per_proc(iproc))
1878 1964 : IF (para_env%mepos == (iproc - 1)) THEN
1879 68979 : INDEX(:) = 0
1880 : counter = 0
1881 4742 : DO i = 1, nhc_len
1882 257321 : DO j = 1, num_nhc
1883 252579 : counter = counter + 1
1884 252579 : work(counter) = nhc%nvt(i, j)%mass
1885 256339 : INDEX(j) = map_info%index(j)
1886 : END DO
1887 : END DO
1888 : ELSE
1889 253561 : work(:) = 0.0_dp
1890 : END IF
1891 1012280 : CALL para_env%bcast(work, iproc - 1)
1892 273952 : CALL para_env%bcast(index, iproc - 1)
1893 1964 : counter = 0
1894 10466 : DO i = 1, nhc_len
1895 514642 : DO j = 1, nhc_per_proc(iproc)
1896 505158 : counter = counter + 1
1897 512678 : mnhc((INDEX(j) - 1)*nhc_len + i) = work(counter)
1898 : END DO
1899 : END DO
1900 : END DO
1901 :
1902 982 : DEALLOCATE (work)
1903 982 : DEALLOCATE (index)
1904 982 : DEALLOCATE (nhc_per_proc)
1905 :
1906 982 : END SUBROUTINE collect_nose_restart_info
1907 :
1908 : ! **************************************************************************************************
1909 : !> \brief routine to dump NEB coordinates and velocities section.. fast implementation
1910 : !> \param coord_section ...
1911 : !> \param array ...
1912 : !> \param narray ...
1913 : !> \param nsize ...
1914 : !> \param nfield ...
1915 : !> \param particle_set ...
1916 : !> \param conv_factor ...
1917 : !> \par History
1918 : !> 12.2006 created [teo]
1919 : !> \author Teodoro Laino
1920 : ! **************************************************************************************************
1921 7148 : SUBROUTINE section_neb_coord_val_set(coord_section, array, narray, nsize, nfield, &
1922 : particle_set, conv_factor)
1923 :
1924 : TYPE(section_vals_type), POINTER :: coord_section
1925 : REAL(KIND=dp), DIMENSION(*) :: array
1926 : INTEGER, INTENT(IN) :: narray, nsize, nfield
1927 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1928 : REAL(KIND=dp) :: conv_factor
1929 :
1930 : INTEGER :: ik, irk, Nlist
1931 7148 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_c
1932 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
1933 : TYPE(section_type), POINTER :: section
1934 : TYPE(val_type), POINTER :: my_val, old_val
1935 :
1936 7148 : NULLIFY (my_val, old_val, section, vals)
1937 0 : CPASSERT(ASSOCIATED(coord_section))
1938 7148 : CPASSERT(coord_section%ref_count > 0)
1939 7148 : section => coord_section%section
1940 7148 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
1941 7148 : IF (ik == -2) &
1942 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
1943 0 : "_DEFAULT_KEYWORD_")
1944 364 : DO
1945 7512 : IF (SIZE(coord_section%values, 2) == 1) EXIT
1946 364 : CALL section_vals_add_values(coord_section)
1947 : END DO
1948 7148 : vals => coord_section%values(ik, 1)%list
1949 7148 : Nlist = 0
1950 7148 : IF (ASSOCIATED(vals)) THEN
1951 6784 : Nlist = cp_sll_val_get_length(vals)
1952 : END IF
1953 270192 : DO irk = 1, nsize/nfield
1954 789132 : ALLOCATE (my_c(nfield))
1955 263044 : IF (nfield == 3) THEN
1956 1051696 : my_c(1:3) = get_particle_pos_or_vel(irk, particle_set, array(1:narray))
1957 1051696 : my_c(1:3) = my_c(1:3)*conv_factor
1958 : ELSE
1959 120 : my_c(1) = array(irk)
1960 : END IF
1961 263044 : CALL val_create(my_val, r_vals_ptr=my_c)
1962 :
1963 263044 : IF (Nlist /= 0) THEN
1964 241140 : IF (irk == 1) THEN
1965 6784 : new_pos => vals
1966 : ELSE
1967 234356 : new_pos => new_pos%rest
1968 : END IF
1969 241140 : old_val => new_pos%first_el
1970 241140 : CALL val_release(old_val)
1971 241140 : new_pos%first_el => my_val
1972 : ELSE
1973 21904 : IF (irk == 1) THEN
1974 364 : NULLIFY (new_pos)
1975 364 : CALL cp_sll_val_create(new_pos, first_el=my_val)
1976 364 : vals => new_pos
1977 : ELSE
1978 21540 : NULLIFY (new_pos%rest)
1979 21540 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
1980 21540 : new_pos => new_pos%rest
1981 : END IF
1982 : END IF
1983 270192 : NULLIFY (my_val)
1984 : END DO
1985 :
1986 7148 : coord_section%values(ik, 1)%list => vals
1987 :
1988 7148 : END SUBROUTINE section_neb_coord_val_set
1989 :
1990 : ! **************************************************************************************************
1991 : !> \brief Set the nose structure like restart
1992 : !> \param work_section ...
1993 : !> \param eta ...
1994 : !> \param veta ...
1995 : !> \param fnhc ...
1996 : !> \param mnhc ...
1997 : !> \par History
1998 : !> 01.2006 created [teo]
1999 : !> \author Teodoro Laino
2000 : ! **************************************************************************************************
2001 1430 : SUBROUTINE set_template_restart(work_section, eta, veta, fnhc, mnhc)
2002 :
2003 : TYPE(section_vals_type), POINTER :: work_section
2004 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eta, veta, fnhc, mnhc
2005 :
2006 : TYPE(section_vals_type), POINTER :: coord, force, mass, velocity
2007 :
2008 1430 : NULLIFY (coord, force, velocity, mass)
2009 1430 : IF (PRESENT(eta)) THEN
2010 1034 : IF (SIZE(eta) > 0) THEN
2011 1034 : coord => section_vals_get_subs_vals(work_section, "COORD")
2012 1034 : CALL section_vals_val_set(coord, "_DEFAULT_KEYWORD_", r_vals_ptr=eta)
2013 : ELSE
2014 0 : DEALLOCATE (eta)
2015 : END IF
2016 : END IF
2017 1430 : IF (PRESENT(veta)) THEN
2018 1430 : IF (SIZE(veta) > 0) THEN
2019 1430 : velocity => section_vals_get_subs_vals(work_section, "VELOCITY")
2020 1430 : CALL section_vals_val_set(velocity, "_DEFAULT_KEYWORD_", r_vals_ptr=veta)
2021 : ELSE
2022 0 : DEALLOCATE (veta)
2023 : END IF
2024 : END IF
2025 1430 : IF (PRESENT(fnhc)) THEN
2026 1034 : IF (SIZE(fnhc) > 0) THEN
2027 1034 : force => section_vals_get_subs_vals(work_section, "FORCE")
2028 1034 : CALL section_vals_val_set(force, "_DEFAULT_KEYWORD_", r_vals_ptr=fnhc)
2029 : ELSE
2030 0 : DEALLOCATE (fnhc)
2031 : END IF
2032 : END IF
2033 1430 : IF (PRESENT(mnhc)) THEN
2034 1430 : IF (SIZE(mnhc) > 0) THEN
2035 1430 : mass => section_vals_get_subs_vals(work_section, "MASS")
2036 1430 : CALL section_vals_val_set(mass, "_DEFAULT_KEYWORD_", r_vals_ptr=mnhc)
2037 : ELSE
2038 0 : DEALLOCATE (mnhc)
2039 : END IF
2040 : END IF
2041 :
2042 1430 : END SUBROUTINE set_template_restart
2043 :
2044 : ! **************************************************************************************************
2045 : !> \brief routine to dump hills information during metadynamics run
2046 : !> \param ss_section ...
2047 : !> \param meta_env ...
2048 : !> \par History
2049 : !> 02.2006 created [teo]
2050 : !> \author Teodoro Laino
2051 : ! **************************************************************************************************
2052 782 : SUBROUTINE meta_hills_val_set_ss(ss_section, meta_env)
2053 :
2054 : TYPE(section_vals_type), POINTER :: ss_section
2055 : TYPE(meta_env_type), POINTER :: meta_env
2056 :
2057 : INTEGER :: ik, irk, lsize, Nlist
2058 782 : REAL(KIND=dp), DIMENSION(:), POINTER :: ss_val
2059 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2060 : TYPE(section_type), POINTER :: section
2061 : TYPE(val_type), POINTER :: my_val, old_val
2062 :
2063 782 : NULLIFY (my_val, old_val, section, vals)
2064 0 : CPASSERT(ASSOCIATED(ss_section))
2065 782 : CPASSERT(ss_section%ref_count > 0)
2066 782 : section => ss_section%section
2067 782 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2068 782 : IF (ik == -2) &
2069 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2070 0 : "_DEFAULT_KEYWORD_")
2071 98 : DO
2072 880 : IF (SIZE(ss_section%values, 2) == 1) EXIT
2073 98 : CALL section_vals_add_values(ss_section)
2074 : END DO
2075 782 : vals => ss_section%values(ik, 1)%list
2076 782 : Nlist = 0
2077 782 : IF (ASSOCIATED(vals)) THEN
2078 684 : Nlist = cp_sll_val_get_length(vals)
2079 : END IF
2080 782 : lsize = SIZE(meta_env%hills_env%ss_history, 1)
2081 12916 : DO irk = 1, meta_env%hills_env%n_hills
2082 36402 : ALLOCATE (ss_val(lsize))
2083 : ! Always stored in A.U.
2084 49176 : ss_val = meta_env%hills_env%ss_history(:, irk)
2085 12134 : CALL val_create(my_val, r_vals_ptr=ss_val)
2086 :
2087 12134 : IF (irk <= Nlist) THEN
2088 10980 : IF (irk == 1) THEN
2089 684 : new_pos => vals
2090 : ELSE
2091 10296 : new_pos => new_pos%rest
2092 : END IF
2093 10980 : old_val => new_pos%first_el
2094 10980 : CALL val_release(old_val)
2095 10980 : new_pos%first_el => my_val
2096 : ELSE
2097 1154 : IF (irk == 1) THEN
2098 98 : NULLIFY (new_pos)
2099 98 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2100 98 : vals => new_pos
2101 : ELSE
2102 1056 : NULLIFY (new_pos%rest)
2103 1056 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2104 1056 : new_pos => new_pos%rest
2105 : END IF
2106 : END IF
2107 12916 : NULLIFY (my_val)
2108 : END DO
2109 :
2110 782 : ss_section%values(ik, 1)%list => vals
2111 :
2112 782 : END SUBROUTINE meta_hills_val_set_ss
2113 :
2114 : ! **************************************************************************************************
2115 : !> \brief routine to dump hills information during metadynamics run
2116 : !> \param ds_section ...
2117 : !> \param meta_env ...
2118 : !> \par History
2119 : !> 02.2006 created [teo]
2120 : !> \author Teodoro Laino
2121 : ! **************************************************************************************************
2122 782 : SUBROUTINE meta_hills_val_set_ds(ds_section, meta_env)
2123 :
2124 : TYPE(section_vals_type), POINTER :: ds_section
2125 : TYPE(meta_env_type), POINTER :: meta_env
2126 :
2127 : INTEGER :: ik, irk, lsize, Nlist
2128 782 : REAL(KIND=dp), DIMENSION(:), POINTER :: ds_val
2129 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2130 : TYPE(section_type), POINTER :: section
2131 : TYPE(val_type), POINTER :: my_val, old_val
2132 :
2133 782 : NULLIFY (my_val, old_val, section, vals)
2134 0 : CPASSERT(ASSOCIATED(ds_section))
2135 782 : CPASSERT(ds_section%ref_count > 0)
2136 782 : section => ds_section%section
2137 782 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2138 782 : IF (ik == -2) &
2139 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2140 0 : "_DEFAULT_KEYWORD_")
2141 98 : DO
2142 880 : IF (SIZE(ds_section%values, 2) == 1) EXIT
2143 98 : CALL section_vals_add_values(ds_section)
2144 : END DO
2145 782 : vals => ds_section%values(ik, 1)%list
2146 782 : Nlist = 0
2147 782 : IF (ASSOCIATED(vals)) THEN
2148 684 : Nlist = cp_sll_val_get_length(vals)
2149 : END IF
2150 782 : lsize = SIZE(meta_env%hills_env%delta_s_history, 1)
2151 12916 : DO irk = 1, meta_env%hills_env%n_hills
2152 36402 : ALLOCATE (ds_val(lsize))
2153 : ! Always stored in A.U.
2154 49176 : ds_val = meta_env%hills_env%delta_s_history(:, irk)
2155 12134 : CALL val_create(my_val, r_vals_ptr=ds_val)
2156 :
2157 12134 : IF (irk <= Nlist) THEN
2158 10980 : IF (irk == 1) THEN
2159 684 : new_pos => vals
2160 : ELSE
2161 10296 : new_pos => new_pos%rest
2162 : END IF
2163 10980 : old_val => new_pos%first_el
2164 10980 : CALL val_release(old_val)
2165 10980 : new_pos%first_el => my_val
2166 : ELSE
2167 1154 : IF (irk == 1) THEN
2168 98 : NULLIFY (new_pos)
2169 98 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2170 98 : vals => new_pos
2171 : ELSE
2172 1056 : NULLIFY (new_pos%rest)
2173 1056 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2174 1056 : new_pos => new_pos%rest
2175 : END IF
2176 : END IF
2177 12916 : NULLIFY (my_val)
2178 : END DO
2179 :
2180 782 : ds_section%values(ik, 1)%list => vals
2181 :
2182 782 : END SUBROUTINE meta_hills_val_set_ds
2183 :
2184 : ! **************************************************************************************************
2185 : !> \brief routine to dump hills information during metadynamics run
2186 : !> \param ww_section ...
2187 : !> \param meta_env ...
2188 : !> \par History
2189 : !> 02.2006 created [teo]
2190 : !> \author Teodoro Laino
2191 : ! **************************************************************************************************
2192 782 : SUBROUTINE meta_hills_val_set_ww(ww_section, meta_env)
2193 :
2194 : TYPE(section_vals_type), POINTER :: ww_section
2195 : TYPE(meta_env_type), POINTER :: meta_env
2196 :
2197 : INTEGER :: ik, irk, lsize, Nlist
2198 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2199 : TYPE(section_type), POINTER :: section
2200 : TYPE(val_type), POINTER :: my_val, old_val
2201 :
2202 782 : NULLIFY (my_val, old_val, section, vals)
2203 782 : CPASSERT(ASSOCIATED(ww_section))
2204 782 : CPASSERT(ww_section%ref_count > 0)
2205 782 : section => ww_section%section
2206 782 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2207 782 : IF (ik == -2) &
2208 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2209 0 : "_DEFAULT_KEYWORD_")
2210 98 : DO
2211 880 : IF (SIZE(ww_section%values, 2) == 1) EXIT
2212 98 : CALL section_vals_add_values(ww_section)
2213 : END DO
2214 782 : vals => ww_section%values(ik, 1)%list
2215 782 : Nlist = 0
2216 782 : IF (ASSOCIATED(vals)) THEN
2217 684 : Nlist = cp_sll_val_get_length(vals)
2218 : END IF
2219 782 : lsize = meta_env%hills_env%n_hills
2220 12916 : DO irk = 1, lsize
2221 12134 : CALL val_create(my_val, r_val=meta_env%hills_env%ww_history(irk))
2222 :
2223 12134 : IF (irk <= Nlist) THEN
2224 10980 : IF (irk == 1) THEN
2225 684 : new_pos => vals
2226 : ELSE
2227 10296 : new_pos => new_pos%rest
2228 : END IF
2229 10980 : old_val => new_pos%first_el
2230 10980 : CALL val_release(old_val)
2231 10980 : new_pos%first_el => my_val
2232 : ELSE
2233 1154 : IF (irk == 1) THEN
2234 98 : NULLIFY (new_pos)
2235 98 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2236 98 : vals => new_pos
2237 : ELSE
2238 1056 : NULLIFY (new_pos%rest)
2239 1056 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2240 1056 : new_pos => new_pos%rest
2241 : END IF
2242 : END IF
2243 12916 : NULLIFY (my_val)
2244 : END DO
2245 :
2246 782 : ww_section%values(ik, 1)%list => vals
2247 :
2248 782 : END SUBROUTINE meta_hills_val_set_ww
2249 :
2250 : ! **************************************************************************************************
2251 : !> \brief routine to dump hills information during metadynamics run
2252 : !> \param invdt_section ...
2253 : !> \param meta_env ...
2254 : !> \par History
2255 : !> 12.2009 created [seb]
2256 : !> \author SC
2257 : ! **************************************************************************************************
2258 2 : SUBROUTINE meta_hills_val_set_dt(invdt_section, meta_env)
2259 :
2260 : TYPE(section_vals_type), POINTER :: invdt_section
2261 : TYPE(meta_env_type), POINTER :: meta_env
2262 :
2263 : INTEGER :: ik, irk, lsize, Nlist
2264 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
2265 : TYPE(section_type), POINTER :: section
2266 : TYPE(val_type), POINTER :: my_val, old_val
2267 :
2268 2 : NULLIFY (my_val, old_val, section, vals)
2269 2 : CPASSERT(ASSOCIATED(invdt_section))
2270 2 : CPASSERT(invdt_section%ref_count > 0)
2271 2 : section => invdt_section%section
2272 2 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
2273 2 : IF (ik == -2) &
2274 : CALL cp_abort(__LOCATION__, "section "//TRIM(section%name)//" does not contain keyword "// &
2275 0 : "_DEFAULT_KEYWORD_")
2276 2 : DO
2277 4 : IF (SIZE(invdt_section%values, 2) == 1) EXIT
2278 2 : CALL section_vals_add_values(invdt_section)
2279 : END DO
2280 2 : vals => invdt_section%values(ik, 1)%list
2281 2 : Nlist = 0
2282 2 : IF (ASSOCIATED(vals)) THEN
2283 0 : Nlist = cp_sll_val_get_length(vals)
2284 : END IF
2285 2 : lsize = meta_env%hills_env%n_hills
2286 6 : DO irk = 1, lsize
2287 4 : CALL val_create(my_val, r_val=meta_env%hills_env%invdt_history(irk))
2288 :
2289 4 : IF (irk <= Nlist) THEN
2290 0 : IF (irk == 1) THEN
2291 0 : new_pos => vals
2292 : ELSE
2293 0 : new_pos => new_pos%rest
2294 : END IF
2295 0 : old_val => new_pos%first_el
2296 0 : CALL val_release(old_val)
2297 0 : new_pos%first_el => my_val
2298 : ELSE
2299 4 : IF (irk == 1) THEN
2300 2 : NULLIFY (new_pos)
2301 2 : CALL cp_sll_val_create(new_pos, first_el=my_val)
2302 2 : vals => new_pos
2303 : ELSE
2304 2 : NULLIFY (new_pos%rest)
2305 2 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
2306 2 : new_pos => new_pos%rest
2307 : END IF
2308 : END IF
2309 6 : NULLIFY (my_val)
2310 : END DO
2311 2 : invdt_section%values(ik, 1)%list => vals
2312 2 : END SUBROUTINE meta_hills_val_set_dt
2313 :
2314 : ! **************************************************************************************************
2315 : !> \brief Write all input sections scaling in size with the number of atoms
2316 : !> in the system to an external file in binary format
2317 : !> \param output_unit binary file to write to
2318 : !> \param log_unit unit for logging debug information
2319 : !> \param root_section ...
2320 : !> \param md_env ...
2321 : !> \param force_env ...
2322 : !> \par History
2323 : !> - Creation (10.02.2011,MK)
2324 : !> \author Matthias Krack (MK)
2325 : !> \version 1.0
2326 : ! **************************************************************************************************
2327 272 : SUBROUTINE write_binary_restart(output_unit, log_unit, root_section, md_env, force_env)
2328 :
2329 : INTEGER, INTENT(IN) :: output_unit, log_unit
2330 : TYPE(section_vals_type), POINTER :: root_section
2331 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
2332 : TYPE(force_env_type), OPTIONAL, POINTER :: force_env
2333 :
2334 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_binary_restart'
2335 :
2336 : CHARACTER(LEN=default_path_length) :: binary_restart_file_name
2337 : CHARACTER(LEN=default_string_length) :: section_label
2338 : INTEGER :: handle, iatom, icore, ikind, imolecule, ishell, istat, n_char_size, n_dp_size, &
2339 : n_int_size, natom, natomkind, ncore, nhc_size, nmolecule, nmoleculekind, nshell, &
2340 : print_level, run_type
2341 272 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf, imol
2342 : LOGICAL :: print_info, write_velocities
2343 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rbuf
2344 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2345 : TYPE(cp_subsys_type), POINTER :: subsys
2346 : TYPE(force_env_type), POINTER :: my_force_env
2347 : TYPE(lnhc_parameters_type), POINTER :: nhc
2348 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2349 : TYPE(molecule_list_type), POINTER :: molecules
2350 : TYPE(mp_para_env_type), POINTER :: para_env
2351 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
2352 : shell_particles
2353 : TYPE(thermostat_type), POINTER :: thermostat_part, thermostat_shell
2354 :
2355 272 : CALL timeset(routineN, handle)
2356 :
2357 272 : NULLIFY (atomic_kinds)
2358 272 : NULLIFY (core_particles)
2359 272 : NULLIFY (molecule_kinds)
2360 272 : NULLIFY (molecules)
2361 272 : NULLIFY (my_force_env)
2362 272 : NULLIFY (para_env)
2363 272 : NULLIFY (particles)
2364 272 : NULLIFY (shell_particles)
2365 272 : NULLIFY (subsys)
2366 272 : NULLIFY (thermostat_part)
2367 272 : NULLIFY (thermostat_shell)
2368 :
2369 272 : IF (PRESENT(md_env)) THEN
2370 : CALL get_md_env(md_env=md_env, &
2371 : force_env=my_force_env, &
2372 : thermostat_part=thermostat_part, &
2373 272 : thermostat_shell=thermostat_shell)
2374 0 : ELSE IF (PRESENT(force_env)) THEN
2375 0 : my_force_env => force_env
2376 : END IF
2377 :
2378 272 : IF (.NOT. ASSOCIATED(my_force_env)) THEN
2379 0 : CALL timestop(handle)
2380 0 : RETURN
2381 : END IF
2382 :
2383 272 : CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", i_val=print_level)
2384 :
2385 272 : IF (print_level > 1) THEN
2386 272 : print_info = .TRUE.
2387 : ELSE
2388 0 : print_info = .FALSE.
2389 : END IF
2390 :
2391 272 : CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=run_type)
2392 : write_velocities = ((run_type == mol_dyn_run) .OR. &
2393 : (run_type == mon_car_run) .OR. &
2394 272 : (run_type == pint_run))
2395 :
2396 : CALL force_env_get(force_env=my_force_env, &
2397 : para_env=para_env, &
2398 272 : subsys=subsys)
2399 : CALL cp_subsys_get(subsys, &
2400 : atomic_kinds=atomic_kinds, &
2401 : particles=particles, &
2402 : natom=natom, &
2403 : core_particles=core_particles, &
2404 : ncore=ncore, &
2405 : shell_particles=shell_particles, &
2406 : nshell=nshell, &
2407 : molecule_kinds=molecule_kinds, &
2408 272 : molecules=molecules)
2409 :
2410 272 : natomkind = atomic_kinds%n_els
2411 272 : IF (ASSOCIATED(molecule_kinds)) THEN
2412 272 : nmoleculekind = molecule_kinds%n_els
2413 : ELSE
2414 0 : nmoleculekind = 0
2415 : END IF
2416 :
2417 272 : IF (ASSOCIATED(molecules)) THEN
2418 272 : nmolecule = molecules%n_els
2419 : ELSE
2420 0 : nmolecule = 0
2421 : END IF
2422 :
2423 272 : n_char_size = 0 ! init
2424 272 : n_int_size = 0 ! init
2425 272 : n_dp_size = 0 ! init
2426 :
2427 272 : IF (output_unit > 0) THEN ! only ionode
2428 :
2429 136 : IF (print_info) THEN
2430 136 : INQUIRE (UNIT=output_unit, NAME=binary_restart_file_name, IOSTAT=istat)
2431 136 : IF (istat /= 0) THEN
2432 : CALL cp_abort(__LOCATION__, &
2433 : "An error occurred inquiring logical unit <"// &
2434 : TRIM(ADJUSTL(cp_to_string(output_unit)))// &
2435 0 : "> which should be linked to the binary restart file")
2436 : END IF
2437 136 : IF (log_unit > 0) THEN
2438 : WRITE (UNIT=log_unit, FMT="(T2,A,/,/,(T3,A,T71,I10))") &
2439 136 : "Writing binary restart file "//TRIM(ADJUSTL(binary_restart_file_name)), &
2440 136 : "Number of atomic kinds:", natomkind, &
2441 136 : "Number of atoms:", natom, &
2442 136 : "Number of cores (only core-shell model):", ncore, &
2443 136 : "Number of shells (only core-shell model):", nshell, &
2444 136 : "Number of molecule kinds:", nmoleculekind, &
2445 272 : "Number of molecules", nmolecule
2446 : END IF
2447 :
2448 136 : n_int_size = n_int_size + 6
2449 : END IF
2450 :
2451 : WRITE (UNIT=output_unit, IOSTAT=istat) &
2452 136 : natomkind, natom, ncore, nshell, nmoleculekind, nmolecule
2453 136 : IF (istat /= 0) THEN
2454 : CALL stop_write("natomkind,natom,ncore,nshell,nmoleculekind,nmolecule "// &
2455 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2456 0 : output_unit)
2457 : END IF
2458 :
2459 : ! Write atomic kind names
2460 408 : DO ikind = 1, natomkind
2461 272 : WRITE (UNIT=output_unit, IOSTAT=istat) atomic_kinds%els(ikind)%name
2462 272 : IF (istat /= 0) CALL stop_write("atomic_kinds%els(ikind)%name "// &
2463 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2464 0 : output_unit)
2465 408 : n_char_size = n_char_size + LEN(atomic_kinds%els(ikind)%name)
2466 : END DO
2467 :
2468 : ! Write atomic kind numbers of all atoms
2469 408 : ALLOCATE (ibuf(natom))
2470 13192 : DO iatom = 1, natom
2471 13192 : ibuf(iatom) = particles%els(iatom)%atomic_kind%kind_number
2472 : END DO
2473 136 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:natom)
2474 136 : IF (istat /= 0) CALL stop_write("ibuf(1:natom) -> atomic kind numbers "// &
2475 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2476 0 : output_unit)
2477 136 : n_int_size = n_int_size + natom
2478 : ! Write atomic coordinates
2479 408 : ALLOCATE (rbuf(3, natom))
2480 13192 : DO iatom = 1, natom
2481 52360 : rbuf(1:3, iatom) = particles%els(iatom)%r(1:3)
2482 : END DO
2483 136 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:natom)
2484 136 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:natom) -> atomic coordinates "// &
2485 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2486 0 : output_unit)
2487 136 : n_dp_size = n_dp_size + 3*natom
2488 136 : DEALLOCATE (rbuf)
2489 :
2490 : ! Write molecule information if available
2491 136 : IF (nmolecule > 0) THEN
2492 : ! Write molecule kind names
2493 272 : DO ikind = 1, nmoleculekind
2494 136 : WRITE (UNIT=output_unit, IOSTAT=istat) molecule_kinds%els(ikind)%name
2495 136 : IF (istat /= 0) CALL stop_write("molecule_kinds%els(ikind)%name "// &
2496 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2497 0 : output_unit)
2498 272 : n_char_size = n_char_size + LEN(molecule_kinds%els(ikind)%name)
2499 : END DO
2500 : ! Write molecule (kind) index numbers for all atoms
2501 13192 : ibuf(:) = 0
2502 408 : ALLOCATE (imol(natom))
2503 13192 : imol(:) = 0
2504 1224 : DO imolecule = 1, nmolecule
2505 1088 : ikind = molecules%els(imolecule)%molecule_kind%kind_number
2506 14144 : DO iatom = molecules%els(imolecule)%first_atom, &
2507 1224 : molecules%els(imolecule)%last_atom
2508 13056 : ibuf(iatom) = ikind
2509 14144 : imol(iatom) = imolecule
2510 : END DO
2511 : END DO
2512 : ! Write molecule kind index number for each atom
2513 136 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:natom)
2514 136 : IF (istat /= 0) CALL stop_write("ibuf(1:natom) -> molecule kind index numbers "// &
2515 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2516 0 : output_unit)
2517 136 : n_int_size = n_int_size + natom
2518 : ! Write molecule index number for each atom
2519 136 : WRITE (UNIT=output_unit, IOSTAT=istat) imol(1:natom)
2520 136 : IF (istat /= 0) CALL stop_write("imol(1:natom) -> molecule index numbers "// &
2521 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2522 0 : output_unit)
2523 136 : n_int_size = n_int_size + natom
2524 136 : DEALLOCATE (imol)
2525 : END IF ! molecules
2526 :
2527 136 : DEALLOCATE (ibuf)
2528 :
2529 : ! Core-shell model only
2530 136 : section_label = "SHELL COORDINATES"
2531 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nshell
2532 136 : IF (istat /= 0) CALL stop_write("section_label, nshell "// &
2533 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2534 0 : output_unit)
2535 136 : n_char_size = n_char_size + LEN(section_label)
2536 136 : n_int_size = n_int_size + 1
2537 136 : IF (nshell > 0) THEN
2538 : ! Write shell coordinates
2539 168 : ALLOCATE (rbuf(3, nshell))
2540 5432 : DO ishell = 1, nshell
2541 21560 : rbuf(1:3, ishell) = shell_particles%els(ishell)%r(1:3)
2542 : END DO
2543 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:nshell)
2544 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:nshell) -> shell coordinates "// &
2545 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2546 0 : output_unit)
2547 56 : n_dp_size = n_dp_size + 3*nshell
2548 56 : DEALLOCATE (rbuf)
2549 : ! Write atomic indices, i.e. number of the atom the shell belongs to
2550 168 : ALLOCATE (ibuf(nshell))
2551 5432 : DO ishell = 1, nshell
2552 5432 : ibuf(ishell) = shell_particles%els(ishell)%atom_index
2553 : END DO
2554 56 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:nshell)
2555 56 : IF (istat /= 0) CALL stop_write("ibuf(1:nshell) -> atomic indices "// &
2556 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2557 0 : output_unit)
2558 56 : n_int_size = n_int_size + nshell
2559 56 : DEALLOCATE (ibuf)
2560 : END IF
2561 :
2562 136 : section_label = "CORE COORDINATES"
2563 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, ncore
2564 136 : IF (istat /= 0) CALL stop_write("section_label, ncore "// &
2565 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2566 0 : output_unit)
2567 136 : n_char_size = n_char_size + LEN(section_label)
2568 136 : n_int_size = n_int_size + 1
2569 136 : IF (ncore > 0) THEN
2570 : ! Write core coordinates
2571 168 : ALLOCATE (rbuf(3, ncore))
2572 5432 : DO icore = 1, ncore
2573 21560 : rbuf(1:3, icore) = core_particles%els(icore)%r(1:3)
2574 : END DO
2575 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:ncore)
2576 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:ncore) -> core coordinates "// &
2577 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2578 0 : output_unit)
2579 56 : n_dp_size = n_dp_size + 3*ncore
2580 56 : DEALLOCATE (rbuf)
2581 : ! Write atomic indices, i.e. number of the atom the core belongs to
2582 168 : ALLOCATE (ibuf(ncore))
2583 5432 : DO icore = 1, ncore
2584 5432 : ibuf(icore) = core_particles%els(icore)%atom_index
2585 : END DO
2586 56 : WRITE (UNIT=output_unit, IOSTAT=istat) ibuf(1:ncore)
2587 56 : IF (istat /= 0) CALL stop_write("ibuf(1:ncore) -> atomic indices "// &
2588 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2589 0 : output_unit)
2590 56 : n_int_size = n_int_size + ncore
2591 56 : DEALLOCATE (ibuf)
2592 : END IF
2593 : END IF ! ionode only
2594 :
2595 : ! Thermostat information
2596 :
2597 : ! Particle thermostats
2598 272 : section_label = "PARTICLE THERMOSTATS"
2599 272 : IF (ASSOCIATED(thermostat_part)) THEN
2600 : ! Nose-Hoover thermostats
2601 176 : IF (thermostat_part%type_of_thermostat == do_thermo_nose) THEN
2602 176 : nhc => thermostat_part%nhc
2603 : CALL write_binary_thermostats_nose(nhc, output_unit, log_unit, section_label, &
2604 : n_char_size, n_dp_size, n_int_size, &
2605 176 : print_info, para_env)
2606 : END IF
2607 : ELSE
2608 96 : nhc_size = 0
2609 96 : IF (output_unit > 0) THEN
2610 48 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nhc_size
2611 48 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", nhc_size "// &
2612 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2613 0 : output_unit)
2614 : END IF
2615 96 : n_char_size = n_char_size + LEN(section_label)
2616 96 : n_int_size = n_int_size + 1
2617 96 : IF (output_unit > 0 .AND. log_unit > 0) THEN ! only ionode
2618 48 : IF (print_info) THEN
2619 : WRITE (UNIT=log_unit, FMT="(T3,A,T71,I10)") &
2620 48 : "NHC size ("//TRIM(ADJUSTL(section_label))//")", nhc_size
2621 : END IF
2622 : END IF
2623 : END IF
2624 :
2625 : ! Shell thermostats (only for core-shell models)
2626 272 : section_label = "SHELL THERMOSTATS"
2627 272 : IF (ASSOCIATED(thermostat_shell)) THEN
2628 : ! Nose-Hoover thermostats
2629 24 : IF (thermostat_shell%type_of_thermostat == do_thermo_nose) THEN
2630 24 : nhc => thermostat_shell%nhc
2631 : CALL write_binary_thermostats_nose(nhc, output_unit, log_unit, section_label, &
2632 : n_char_size, n_dp_size, n_int_size, &
2633 24 : print_info, para_env)
2634 : END IF
2635 : ELSE
2636 248 : nhc_size = 0
2637 248 : IF (output_unit > 0) THEN
2638 124 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nhc_size
2639 124 : IF (istat /= 0) CALL stop_write("nhc_size "// &
2640 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2641 0 : output_unit)
2642 : END IF
2643 248 : n_char_size = n_char_size + LEN(section_label)
2644 248 : n_int_size = n_int_size + 1
2645 248 : IF (output_unit > 0 .AND. log_unit > 0) THEN ! only ionode
2646 124 : IF (print_info) THEN
2647 : WRITE (UNIT=log_unit, FMT="(T3,A,T71,I10)") &
2648 124 : "NHC size ("//TRIM(ADJUSTL(section_label))//")", nhc_size
2649 : END IF
2650 : END IF
2651 : END IF
2652 :
2653 : ! Particle velocities
2654 :
2655 272 : IF (output_unit > 0) THEN ! only ionode
2656 : ! Write particle velocities if needed
2657 136 : section_label = "VELOCITIES"
2658 : IF (output_unit > 0) THEN
2659 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, MERGE(natom, 0, write_velocities)
2660 136 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", write_velocities "// &
2661 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2662 0 : output_unit)
2663 : END IF
2664 136 : n_char_size = n_char_size + LEN(section_label)
2665 136 : n_int_size = n_int_size + 1
2666 136 : IF (print_info .AND. log_unit > 0) THEN
2667 : WRITE (UNIT=log_unit, FMT="(T3,A,T78,A3)") &
2668 136 : "Write "//TRIM(ADJUSTL(section_label))//" section", MERGE("YES", " NO", write_velocities)
2669 : END IF
2670 136 : IF (write_velocities) THEN
2671 408 : ALLOCATE (rbuf(3, natom))
2672 : ! Write atomic velocities
2673 13192 : DO iatom = 1, natom
2674 52360 : rbuf(1:3, iatom) = particles%els(iatom)%v(1:3)
2675 : END DO
2676 136 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:natom)
2677 136 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:natom) -> atomic velocities "// &
2678 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2679 0 : output_unit)
2680 136 : n_dp_size = n_dp_size + 3*natom
2681 136 : DEALLOCATE (rbuf)
2682 : END IF
2683 : ! Write shell velocities
2684 136 : section_label = "SHELL VELOCITIES"
2685 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, MERGE(nshell, 0, write_velocities)
2686 136 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", write_velocities "// &
2687 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2688 0 : output_unit)
2689 136 : n_char_size = n_char_size + LEN(section_label)
2690 136 : n_int_size = n_int_size + 1
2691 136 : IF (print_info .AND. log_unit > 0) THEN
2692 : WRITE (UNIT=log_unit, FMT="(T3,A,T78,A3)") &
2693 136 : "Write "//TRIM(ADJUSTL(section_label))//" section", MERGE("YES", " NO", write_velocities)
2694 : END IF
2695 136 : IF (nshell > 0) THEN
2696 56 : IF (write_velocities) THEN
2697 168 : ALLOCATE (rbuf(3, nshell))
2698 5432 : DO ishell = 1, nshell
2699 21560 : rbuf(1:3, ishell) = shell_particles%els(ishell)%v(1:3)
2700 : END DO
2701 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:nshell)
2702 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:nshell) -> shell velocities "// &
2703 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2704 0 : output_unit)
2705 56 : n_dp_size = n_dp_size + 3*nshell
2706 56 : DEALLOCATE (rbuf)
2707 : END IF
2708 : END IF
2709 : ! Write core velocities
2710 136 : section_label = "CORE VELOCITIES"
2711 136 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, MERGE(ncore, 0, write_velocities)
2712 136 : IF (istat /= 0) CALL stop_write(TRIM(section_label)//", write_velocities "// &
2713 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2714 0 : output_unit)
2715 136 : n_char_size = n_char_size + LEN(section_label)
2716 136 : n_int_size = n_int_size + 1
2717 136 : IF (print_info .AND. log_unit > 0) THEN
2718 : WRITE (UNIT=log_unit, FMT="(T3,A,T78,A3)") &
2719 136 : "Write "//TRIM(ADJUSTL(section_label))//" section", MERGE("YES", " NO", write_velocities)
2720 : END IF
2721 136 : IF (ncore > 0) THEN
2722 56 : IF (write_velocities) THEN
2723 168 : ALLOCATE (rbuf(3, ncore))
2724 5432 : DO icore = 1, ncore
2725 21560 : rbuf(1:3, icore) = core_particles%els(icore)%v(1:3)
2726 : END DO
2727 56 : WRITE (UNIT=output_unit, IOSTAT=istat) rbuf(1:3, 1:ncore)
2728 56 : IF (istat /= 0) CALL stop_write("rbuf(1:3,1:ncore) -> core velocities "// &
2729 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2730 0 : output_unit)
2731 56 : n_dp_size = n_dp_size + 3*ncore
2732 56 : DEALLOCATE (rbuf)
2733 : END IF
2734 : END IF
2735 : END IF ! ionode only
2736 :
2737 : ! Optionally, print a small I/O statistics
2738 272 : IF (output_unit > 0) THEN ! only ionode
2739 136 : IF (print_info .AND. log_unit > 0) THEN
2740 : WRITE (UNIT=log_unit, FMT="(/,(T2,I10,1X,I0,A,T68,I10,A))") &
2741 136 : n_char_size, int_size, "-byte characters written", n_char_size*int_size/1024, " KB", &
2742 136 : n_dp_size, dp_size, "-byte floating point numbers written", n_dp_size*dp_size/1024, " KB", &
2743 272 : n_int_size, int_size, "-byte integer numbers written", n_int_size*int_size/1024, " KB"
2744 : WRITE (UNIT=log_unit, FMT="(/,T2,A)") &
2745 136 : "Binary restart file "//TRIM(ADJUSTL(binary_restart_file_name))//" written"
2746 : END IF
2747 : END IF ! ionode only
2748 :
2749 272 : CALL timestop(handle)
2750 :
2751 : END SUBROUTINE write_binary_restart
2752 :
2753 : ! **************************************************************************************************
2754 : !> \brief Write an input section for Nose thermostats to an external file in
2755 : !> binary format
2756 : !> \param nhc ...
2757 : !> \param output_unit binary file to write to
2758 : !> \param log_unit unit for logging debug information
2759 : !> \param section_label ...
2760 : !> \param n_char_size ...
2761 : !> \param n_dp_size ...
2762 : !> \param n_int_size ...
2763 : !> \param print_info ...
2764 : !> \param para_env ...
2765 : !> \par History
2766 : !> - Creation (23.03.2011,MK)
2767 : !> \author Matthias Krack (MK)
2768 : !> \version 1.0
2769 : ! **************************************************************************************************
2770 200 : SUBROUTINE write_binary_thermostats_nose(nhc, output_unit, log_unit, section_label, &
2771 : n_char_size, n_dp_size, n_int_size, &
2772 : print_info, para_env)
2773 :
2774 : TYPE(lnhc_parameters_type), POINTER :: nhc
2775 : INTEGER, INTENT(IN) :: output_unit, log_unit
2776 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section_label
2777 : INTEGER, INTENT(INOUT) :: n_char_size, n_dp_size, n_int_size
2778 : LOGICAL, INTENT(IN) :: print_info
2779 : TYPE(mp_para_env_type), POINTER :: para_env
2780 :
2781 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_binary_thermostats_nose'
2782 :
2783 : INTEGER :: handle, istat, nhc_size
2784 200 : REAL(KIND=dp), DIMENSION(:), POINTER :: eta, fnhc, mnhc, veta
2785 :
2786 200 : CALL timeset(routineN, handle)
2787 :
2788 200 : NULLIFY (eta)
2789 200 : NULLIFY (fnhc)
2790 200 : NULLIFY (mnhc)
2791 200 : NULLIFY (veta)
2792 :
2793 200 : CALL collect_nose_restart_info(nhc, para_env, eta, veta, fnhc, mnhc)
2794 :
2795 200 : nhc_size = SIZE(eta)
2796 :
2797 200 : IF (output_unit > 0) THEN ! only ionode
2798 100 : WRITE (UNIT=output_unit, IOSTAT=istat) section_label, nhc_size
2799 100 : IF (istat /= 0) CALL stop_write("nhc_size "// &
2800 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2801 0 : output_unit)
2802 100 : n_char_size = n_char_size + LEN(section_label)
2803 100 : n_int_size = n_int_size + 1
2804 100 : IF (print_info .AND. log_unit > 0) THEN
2805 : WRITE (UNIT=log_unit, FMT="(T3,A,T71,I10)") &
2806 100 : "NHC size ("//TRIM(ADJUSTL(section_label))//")", nhc_size
2807 : END IF
2808 : ! eta
2809 100 : WRITE (UNIT=output_unit, IOSTAT=istat) eta(1:nhc_size)
2810 100 : IF (istat /= 0) CALL stop_write("eta(1:nhc_size) "// &
2811 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2812 0 : output_unit)
2813 100 : n_dp_size = n_dp_size + nhc_size
2814 : END IF ! ionode only
2815 :
2816 200 : DEALLOCATE (eta)
2817 :
2818 : ! veta
2819 200 : IF (output_unit > 0) THEN ! only ionode
2820 100 : WRITE (UNIT=output_unit, IOSTAT=istat) veta(1:nhc_size)
2821 100 : IF (istat /= 0) CALL stop_write("veta(1:nhc_size) "// &
2822 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2823 0 : output_unit)
2824 100 : n_dp_size = n_dp_size + nhc_size
2825 : END IF ! ionode only
2826 :
2827 200 : DEALLOCATE (veta)
2828 :
2829 : ! mnhc
2830 200 : IF (output_unit > 0) THEN ! only ionode
2831 100 : WRITE (UNIT=output_unit, IOSTAT=istat) mnhc(1:nhc_size)
2832 100 : IF (istat /= 0) CALL stop_write("mnhc(1:nhc_size) "// &
2833 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2834 0 : output_unit)
2835 100 : n_dp_size = n_dp_size + nhc_size
2836 : END IF ! ionode only
2837 :
2838 200 : DEALLOCATE (mnhc)
2839 :
2840 : ! fnhc
2841 200 : IF (output_unit > 0) THEN ! only ionode
2842 100 : WRITE (UNIT=output_unit, IOSTAT=istat) fnhc(1:nhc_size)
2843 100 : IF (istat /= 0) CALL stop_write("fnhc(1:nhc_size) "// &
2844 : "(IOSTAT = "//TRIM(ADJUSTL(cp_to_string(istat)))//")", &
2845 0 : output_unit)
2846 100 : n_dp_size = n_dp_size + nhc_size
2847 : END IF ! ionode only
2848 :
2849 200 : DEALLOCATE (fnhc)
2850 :
2851 200 : CALL timestop(handle)
2852 :
2853 200 : END SUBROUTINE write_binary_thermostats_nose
2854 :
2855 : ! **************************************************************************************************
2856 : !> \brief Print an error message and stop the program execution in case of a
2857 : !> read error.
2858 : !> \param object ...
2859 : !> \param unit_number ...
2860 : !> \par History
2861 : !> - Creation (15.02.2011,MK)
2862 : !> \author Matthias Krack (MK)
2863 : !> \note
2864 : !> object : Name of the data object for which I/O operation failed
2865 : !> unit_number: Logical unit number of the file written to
2866 : ! **************************************************************************************************
2867 0 : SUBROUTINE stop_write(object, unit_number)
2868 :
2869 : CHARACTER(LEN=*), INTENT(IN) :: object
2870 : INTEGER, INTENT(IN) :: unit_number
2871 :
2872 : CHARACTER(LEN=2*default_path_length) :: message
2873 : CHARACTER(LEN=default_path_length) :: file_name
2874 : LOGICAL :: file_exists
2875 :
2876 0 : IF (unit_number >= 0) THEN
2877 0 : INQUIRE (UNIT=unit_number, EXIST=file_exists)
2878 : ELSE
2879 0 : file_exists = .FALSE.
2880 : END IF
2881 0 : IF (file_exists) THEN
2882 0 : INQUIRE (UNIT=unit_number, NAME=file_name)
2883 : WRITE (UNIT=message, FMT="(A)") &
2884 : "An error occurred writing data object <"//TRIM(ADJUSTL(object))// &
2885 0 : "> to file <"//TRIM(ADJUSTL(file_name))//">"
2886 : ELSE
2887 : WRITE (UNIT=message, FMT="(A,I0,A)") &
2888 : "Could not write data object <"//TRIM(ADJUSTL(object))// &
2889 0 : "> to logical unit ", unit_number, ". The I/O unit does not exist."
2890 : END IF
2891 :
2892 0 : CPABORT(message)
2893 :
2894 0 : END SUBROUTINE stop_write
2895 :
2896 : END MODULE input_cp2k_restarts
|
