LCOV - code coverage report
Current view: top level - src/motion - input_cp2k_restarts.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 1278 1352 94.5 %
Date: 2023-03-30 11:55:16 Functions: 20 21 95.2 %

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

Generated by: LCOV version 1.15