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