LCOV - code coverage report
Current view: top level - src/motion - input_cp2k_restarts.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 1243 1352 91.9 %
Date: 2024-04-20 06:29:22 Functions: 19 21 90.5 %

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

Generated by: LCOV version 1.15