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