tamc_run.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Perform a temperature accelarated hybrid monte carlo (TAHMC) run using QUICKSTEP
!> \par History
!>      none
!> \author Alin M Elena
! **************************************************************************************************
MODULE tamc_run
 
   USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
   USE atomic_kind_types,               ONLY: atomic_kind_type
   USE averages_types,                  ONLY: average_quantities_type
   USE barostat_types,                  ONLY: barostat_type,&
                                              create_barostat_type
   USE bibliography,                    ONLY: vandencic2006
   USE cell_types,                      ONLY: cell_type
   USE colvar_methods,                  ONLY: colvar_eval_glob_f
   USE colvar_types,                    ONLY: hbp_colvar_id,&
                                              wc_colvar_id,&
                                              colvar_p_type
   USE constraint_fxd,                  ONLY: fix_atom_control
   USE cp_external_control,             ONLY: external_control
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_get_default_io_unit,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_add_iter_level,&
                                              cp_iterate,&
                                              cp_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr,&
                                              cp_rm_iter_level
   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
                                              cp_subsys_type
   USE cp_units,                        ONLY: cp_unit_from_cp2k
   USE distribution_1d_types,           ONLY: distribution_1d_type
   USE force_env_methods,               ONLY: force_env_calc_energy_force
   USE force_env_types,                 ONLY: force_env_get,&
                                              force_env_type
   USE free_energy_types,               ONLY: fe_env_create,&
                                              free_energy_type
   USE global_types,                    ONLY: global_environment_type
   USE input_constants,                 ONLY: &
        langevin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
        nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, reftraj_ensemble
   USE input_cp2k_check,                ONLY: remove_restart_info
   USE input_cp2k_restarts,             ONLY: write_restart
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_remove_values,&
                                              section_vals_type,&
                                              section_vals_val_get,&
                                              section_vals_val_set
   USE kinds,                           ONLY: dp
   USE machine,                         ONLY: m_walltime
   USE mc_environment_types,            ONLY: get_mc_env,&
                                              mc_env_create,&
                                              mc_env_release,&
                                              mc_environment_type,&
                                              set_mc_env
   USE mc_misc,                         ONLY: mc_averages_create,&
                                              mc_averages_release
   USE mc_move_control,                 ONLY: init_mc_moves,&
                                              mc_moves_release
   USE mc_types,                        ONLY: get_mc_par,&
                                              mc_averages_type,&
                                              mc_ekin_type,&
                                              mc_moves_type,&
                                              mc_simpar_type,&
                                              set_mc_par
   USE md_ener_types,                   ONLY: create_md_ener,&
                                              md_ener_type
   USE md_energies,                     ONLY: initialize_md_ener,&
                                              md_energy
   USE md_environment_types,            ONLY: get_md_env,&
                                              md_env_create,&
                                              md_env_release,&
                                              md_environment_type,&
                                              set_md_env
   USE md_run,                          ONLY: qs_mol_dyn
   USE message_passing,                 ONLY: mp_comm_type,&
                                              mp_para_env_type
   USE metadynamics_types,              ONLY: meta_env_type,&
                                              metavar_type,&
                                              set_meta_env
   USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
   USE molecule_kind_types,             ONLY: molecule_kind_type
   USE molecule_list_types,             ONLY: molecule_list_type
   USE molecule_types,                  ONLY: global_constraint_type,&
                                              molecule_type
   USE parallel_rng_types,              ONLY: uniform,&
                                              rng_stream_type
   USE particle_list_types,             ONLY: particle_list_type
   USE particle_types,                  ONLY: particle_type
   USE physcon,                         ONLY: boltzmann,&
                                              femtoseconds,&
                                              joule,&
                                              kelvin
   USE qmmm_util,                       ONLY: apply_qmmm_walls_reflective
   USE qs_environment_types,            ONLY: get_qs_env
   USE qs_scf_post_gpw,                 ONLY: scf_post_calculation_gpw
   USE reference_manager,               ONLY: cite_reference
   USE reftraj_types,                   ONLY: create_reftraj,&
                                              reftraj_type
   USE reftraj_util,                    ONLY: initialize_reftraj
   USE simpar_methods,                  ONLY: read_md_section
   USE simpar_types,                    ONLY: create_simpar_type,&
                                              release_simpar_type,&
                                              simpar_type
   USE string_utilities,                ONLY: str_comp
   USE thermal_region_types,            ONLY: thermal_regions_type
   USE thermal_region_utils,            ONLY: create_thermal_regions
   USE thermostat_methods,              ONLY: create_thermostats
   USE thermostat_types,                ONLY: thermostats_type
   USE virial_methods,                  ONLY: virial_evaluate
   USE virial_types,                    ONLY: virial_type
   USE wiener_process,                  ONLY: create_wiener_process,&
                                              create_wiener_process_cv
!!!!! monte carlo part
#include "../../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tamc_run'
 
   PUBLIC :: qs_tamc
 
CONTAINS
 
! **************************************************************************************************
!> \brief Driver routine for TAHMC
!> \param force_env ...
!> \param globenv ...
!> \param averages ...
!> \author Alin M Elena
!> \note it computes the forces using QuickStep.
! **************************************************************************************************
   SUBROUTINE qs_tamc(force_env, globenv, averages)
 
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(global_environment_type), POINTER             :: globenv
      TYPE(average_quantities_type), OPTIONAL, POINTER   :: averages
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'qs_tamc'
 
      CHARACTER(LEN=20)                                  :: ensemble
      INTEGER                                            :: handle, i, initialstep, iprint, isos, &
                                                            istep, j, md_stride, nmccycles, &
                                                            output_unit, rand2skip, run_type_id
      INTEGER, POINTER                                   :: itimes
      LOGICAL                                            :: check, explicit, my_rm_restart_info, &
                                                            save_mem, should_stop
      REAL(kind=dp)                                      :: auxrandom, inittime, rval, temp, &
                                                            time_iter_start, time_iter_stop
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: an, fz, xieta, zbuff
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: r
      REAL(kind=dp), POINTER                             :: constant, time, used_time
      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      TYPE(barostat_type), POINTER                       :: barostat
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(cp_subsys_type), POINTER                      :: subsys, subsys_i
      TYPE(distribution_1d_type), POINTER                :: local_particles
      TYPE(free_energy_type), POINTER                    :: fe_env
      TYPE(mc_averages_type), POINTER                    :: mcaverages
      TYPE(mc_environment_type), POINTER                 :: mc_env
      TYPE(mc_moves_type), POINTER                       :: gmoves, moves
      TYPE(mc_simpar_type), POINTER                      :: mc_par
      TYPE(md_ener_type), POINTER                        :: md_ener
      TYPE(md_environment_type), POINTER                 :: md_env
      TYPE(meta_env_type), POINTER                       :: meta_env_saved
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(reftraj_type), POINTER                        :: reftraj
      TYPE(rng_stream_type)                              :: rng_stream_mc
      TYPE(section_vals_type), POINTER :: constraint_section, force_env_section, &
         free_energy_section, fs_section, global_section, mc_section, md_section, motion_section, &
         reftraj_section, subsys_section, work_section
      TYPE(simpar_type), POINTER                         :: simpar
      TYPE(thermal_regions_type), POINTER                :: thermal_regions
      TYPE(thermostats_type), POINTER                    :: thermostats
      TYPE(virial_type), POINTER                         :: virial
 
      initialstep = 0
      inittime = 0.0_dp
 
      CALL timeset(routinen, handle)
      my_rm_restart_info = .true.
      NULLIFY (para_env, fs_section, virial)
      para_env => force_env%para_env
      motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
      md_section => section_vals_get_subs_vals(motion_section, "MD")
 
      ! Real call to MD driver - Low Level
      ALLOCATE (md_env)
      CALL md_env_create(md_env, md_section, para_env, force_env=force_env)
      IF (PRESENT(averages)) CALL set_md_env(md_env, averages=averages)
 
      cpassert(ASSOCIATED(globenv))
      cpassert(ASSOCIATED(force_env))
 
      NULLIFY (particles, cell, simpar, itimes, used_time, subsys, &
               md_ener, thermostats, barostat, reftraj, force_env_section, &
               reftraj_section, work_section, atomic_kinds, &
               local_particles, time, fe_env, free_energy_section, &
               constraint_section, thermal_regions, subsys_i)
      logger => cp_get_default_logger()
      para_env => force_env%para_env
 
      global_section => section_vals_get_subs_vals(force_env%root_section, "GLOBAL")
      free_energy_section => section_vals_get_subs_vals(motion_section, "FREE_ENERGY")
      constraint_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT")
      CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
 
      CALL section_vals_val_get(global_section, "RUN_TYPE", i_val=run_type_id)
 
      CALL create_simpar_type(simpar)
      force_env_section => force_env%force_env_section
      subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
      CALL cp_add_iter_level(logger%iter_info, "MD")
      CALL cp_iterate(logger%iter_info, iter_nr=initialstep)
      ! Read MD section
      CALL read_md_section(simpar, motion_section, md_section)
      ! Setup print_keys
      simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, &
                                                    "CONSTRAINT_INFO", extension=".shakeLog", log_filename=.false.)
      simpar%lagrange_multipliers = cp_print_key_unit_nr(logger, constraint_section, &
                                                         "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.false.)
      simpar%dump_lm = btest(cp_print_key_should_output(logger%iter_info, constraint_section, &
                                                        "LAGRANGE_MULTIPLIERS"), cp_p_file)
 
      ! Create the structure for the md energies
      ALLOCATE (md_ener)
      CALL create_md_ener(md_ener)
      CALL set_md_env(md_env, md_ener=md_ener)
 
      ! If requested setup Thermostats
      CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env, &
                              globenv, global_section)
 
      ! If requested setup Barostat
      CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv)
 
      ! If requested setup different thermal regions
      CALL create_thermal_regions(thermal_regions, md_section, simpar, force_env)
 
      CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions)
 
      IF (simpar%ensemble == reftraj_ensemble) THEN
         reftraj_section => section_vals_get_subs_vals(md_section, "REFTRAJ")
         ALLOCATE (reftraj)
         CALL create_reftraj(reftraj, reftraj_section, para_env)
         CALL set_md_env(md_env, reftraj=reftraj)
      END IF
 
      CALL force_env_get(force_env, subsys=subsys, cell=cell, &
                         force_env_section=force_env_section)
 
      ! Set V0 if needed
      IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
         IF (simpar%v0 == 0._dp) simpar%v0 = cell%deth
      END IF
 
      ! Initialize velocities possibly applying constraints at the zeroth MD step
! ! !     CALL section_vals_val_get(motion_section,"PRINT%RESTART%SPLIT_RESTART_FILE",&
! ! !                               l_val=write_binary_restart_file)
!! let us see if this created all the trouble
!      CALL setup_velocities(force_env,simpar,globenv,md_env,md_section,constraint_section, &
!                            write_binary_restart_file)
 
      ! Setup Free Energy Calculation (if required)
      CALL fe_env_create(fe_env, free_energy_section)
      CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell, &
                      force_env=force_env)
 
      ! Possibly initialize Wiener processes
      IF (simpar%ensemble == langevin_ensemble) CALL create_wiener_process(md_env)
      time_iter_start = m_walltime()
 
      CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, &
                      md_ener=md_ener, t=time, used_time=used_time)
 
      ! Attach the time counter of the meta_env to the one of the MD
      CALL set_meta_env(force_env%meta_env, time=time)
      ! Initialize the md_ener structure
 
      force_env%meta_env%dt = force_env%meta_env%zdt
      CALL initialize_md_ener(md_ener, force_env, simpar)
!     force_env%meta_env%dt=force_env%meta_env%zdt
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MC setup up
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
      NULLIFY (mc_env, mc_par, mcaverages)
 
      CALL section_vals_get(force_env_section, n_repetition=isos)
      cpassert(isos == 1)
! set some values...will use get_globenv if that ever comes around
 
! initialize the random numbers
!       IF (para_env%is_source()) THEN
      rng_stream_mc = rng_stream_type(name="Random numbers for monte carlo acc/rej", &
                                      distribution_type=uniform)
!       ENDIF
!!!!! this should go in a routine hmc_read
 
      NULLIFY (mc_section)
      ALLOCATE (mc_par)
 
      mc_section => section_vals_get_subs_vals(force_env%root_section, &
                                               "MOTION%MC")
      CALL section_vals_val_get(mc_section, "ENSEMBLE", &
                                c_val=ensemble)
      cpassert(str_comp(ensemble, "TRADITIONAL"))
      CALL section_vals_val_get(mc_section, "NSTEP", &
                                i_val=nmccycles)
      cpassert(nmccycles > 0)
      CALL section_vals_val_get(mc_section, "IPRINT", &
                                i_val=iprint)
      CALL section_vals_val_get(mc_section, "RANDOMTOSKIP", i_val=rand2skip)
      cpassert(rand2skip >= 0)
      temp = cp_unit_from_cp2k(simpar%temp_ext, "K")
!
 
      CALL set_mc_par(mc_par, ensemble=ensemble, nstep=nmccycles, iprint=iprint, temperature=temp, &
                      beta=1.0_dp/temp/boltzmann*joule, exp_max_val=0.9_dp*log(huge(0.0_dp)), &
                      exp_min_val=0.9_dp*log(tiny(0.0_dp)), max_val=huge(0.0_dp), min_val=0.0_dp, &
                      source=para_env%source, group=para_env, ionode=para_env%is_source(), rand2skip=rand2skip)
 
      output_unit = cp_logger_get_default_io_unit(logger)
      IF (output_unit > 0) THEN
         WRITE (output_unit, '(a,a)') "HMC| Hybrid Monte Carlo Scheme "
         WRITE (output_unit, '(a,a)') "HMC| Ensemble ", adjustl(ensemble)
         WRITE (output_unit, '(a,i0)') "HMC| MC Cycles ", nmccycles
         WRITE (output_unit, '(a,i0,a)') "HMC| Print every ", iprint, " cycles"
         WRITE (output_unit, '(a,i0)') "HMC| Number of random numbers to skip ", rand2skip
         WRITE (output_unit, '(a,f16.8,a)') "HMC| Temperature ", temp, "K"
      END IF
 
      CALL force_env_get(force_env, subsys=subsys)
 
      CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
                         particles=particles, virial=virial)
 
      DO i = 1, rand2skip
         auxrandom = rng_stream_mc%next()
         DO j = 1, 3*SIZE(particles%els)
            auxrandom = globenv%gaussian_rng_stream%next()
         END DO
      END DO
 
      ALLOCATE (mc_env)
      CALL mc_env_create(mc_env)
      CALL set_mc_env(mc_env, mc_par=mc_par, force_env=force_env)
!!!!!!!end mc setup
 
      ! Check for ensembles requiring the stress tensor - takes into account the possibility for
      ! multiple force_evals
      IF ((simpar%ensemble == npt_i_ensemble) .OR. &
          (simpar%ensemble == npt_ia_ensemble) .OR. &
          (simpar%ensemble == npt_f_ensemble) .OR. &
          (simpar%ensemble == npe_f_ensemble) .OR. &
          (simpar%ensemble == npe_i_ensemble) .OR. &
          (simpar%ensemble == nph_uniaxial_ensemble) .OR. &
          (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
         check = virial%pv_availability
         IF (.NOT. check) &
            CALL cp_abort(__location__, &
                          "Virial evaluation not requested for this run in the input file! "// &
                          "You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR. "// &
                          "Be sure the method you are using can compute the virial!")
         IF (ASSOCIATED(force_env%sub_force_env)) THEN
            DO i = 1, SIZE(force_env%sub_force_env)
               IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN
                  CALL force_env_get(force_env%sub_force_env(i)%force_env, subsys=subsys_i)
                  CALL cp_subsys_get(subsys_i, virial=virial)
                  check = check .AND. virial%pv_availability
               END IF
            END DO
         END IF
         IF (.NOT. check) &
            CALL cp_abort(__location__, &
                          "Virial evaluation not requested for all the force_eval sections present in"// &
                          " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR"// &
                          " in each force_eval section. Be sure the method you are using can compute the virial!")
      END IF
 
      ! Computing Forces at zero MD step
      IF (simpar%ensemble /= reftraj_ensemble) THEN
         CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
         CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time)
         CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant)
         CALL section_vals_val_set(md_section, "STEP_START_VAL", i_val=initialstep)
         CALL section_vals_val_set(md_section, "TIME_START_VAL", r_val=inittime)
         initialstep = itimes
         CALL cp_iterate(logger%iter_info, iter_nr=itimes)
         IF (save_mem) THEN
            work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
            CALL section_vals_remove_values(work_section)
            work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
            CALL section_vals_remove_values(work_section)
            work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
            CALL section_vals_remove_values(work_section)
         END IF
 
!      CALL force_env_calc_energy_force (force_env, calc_force=.TRUE.)
         meta_env_saved => force_env%meta_env
         NULLIFY (force_env%meta_env)
         CALL force_env_calc_energy_force(force_env, calc_force=.false.)
         force_env%meta_env => meta_env_saved
 
         IF (ASSOCIATED(force_env%qs_env)) THEN
!           force_env%qs_env%sim_time=time
!           force_env%qs_env%sim_step=itimes
            force_env%qs_env%sim_time = 0.0_dp
            force_env%qs_env%sim_step = 0
         END IF
         ! Warm-up engines for metadynamics
         IF (ASSOCIATED(force_env%meta_env)) THEN
            IF (force_env%meta_env%langevin) THEN
               CALL create_wiener_process_cv(force_env%meta_env)
               DO j = 1, (rand2skip - 1)/nmccycles
                  DO i = 1, force_env%meta_env%n_colvar
                     auxrandom = force_env%meta_env%rng(i)%next()
                     auxrandom = force_env%meta_env%rng(i)%next()
                  END DO
               END DO
            END IF
!           IF (force_env%meta_env%well_tempered) THEN
!              force_env%meta_env%wttemperature = simpar%temp_ext
!              IF (force_env%meta_env%wtgamma>EPSILON(1._dp)) THEN
!                 dummy=force_env%meta_env%wttemperature*(force_env%meta_env%wtgamma-1._dp)
!                 IF (force_env%meta_env%delta_t>EPSILON(1._dp)) THEN
!                    check=ABS(force_env%meta_env%delta_t-dummy)<1.E+3_dp*EPSILON(1._dp)
!                    IF(.NOT.check) CALL cp_abort(__LOCATION__,&
!                       "Inconsistency between DELTA_T and WTGAMMA (both specified):"//&
!                       " please, verify that DELTA_T=(WTGAMMA-1)*TEMPERATURE")
!                 ELSE
!                    force_env%meta_env%delta_t = dummy
!                 ENDIF
!              ELSE
!                 force_env%meta_env%wtgamma    = 1._dp &
!                    + force_env%meta_env%delta_t/force_env%meta_env%wttemperature
!              ENDIF
!              force_env%meta_env%invdt         = 1._dp/force_env%meta_env%delta_t
!           ENDIF
            CALL tamc_force(force_env)
!           CALL metadyn_write_colvar(force_env)
         END IF
 
         IF (simpar%do_respa) THEN
            CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
                                             calc_force=.true.)
         END IF
 
!        CALL force_env_get( force_env, subsys=subsys)
!
!        CALL cp_subsys_get(subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,&
!             particles=particles)
 
         CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
                              virial, force_env%para_env)
 
         CALL md_energy(md_env, md_ener)
!        CALL md_write_output(md_env) !inits the print env at itimes == 0 also writes trajectories
         md_stride = 1
      ELSE
         CALL get_md_env(md_env, reftraj=reftraj)
         CALL initialize_reftraj(reftraj, reftraj_section, md_env)
         itimes = reftraj%info%first_snapshot - 1
         md_stride = reftraj%info%stride
      END IF
 
      CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
                                        constraint_section, "CONSTRAINT_INFO")
      CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
                                        constraint_section, "LAGRANGE_MULTIPLIERS")
      CALL init_mc_moves(moves)
      CALL init_mc_moves(gmoves)
      ALLOCATE (r(1:3, SIZE(particles%els)))
!       ALLOCATE (r_old(1:3,size(particles%els)))
      CALL mc_averages_create(mcaverages)
      !!!!! some more buffers
      ! Allocate random number for Langevin Thermostat acting on COLVARS
      ALLOCATE (xieta(2*force_env%meta_env%n_colvar))
      xieta(:) = 0.0_dp
      ALLOCATE (an(force_env%meta_env%n_colvar))
      an(:) = 0.0_dp
      ALLOCATE (fz(force_env%meta_env%n_colvar))
      fz(:) = 0.0_dp
      ALLOCATE (zbuff(2*force_env%meta_env%n_colvar))
      zbuff(:) = 0.0_dp
 
      IF (output_unit > 0) THEN
         WRITE (output_unit, '(a)') "HMC|==== Initial average forces"
      END IF
      CALL metadyn_write_colvar_header(force_env)
      moves%hmc%attempts = 0
      moves%hmc%successes = 0
      gmoves%hmc%attempts = 0
      gmoves%hmc%successes = 0
      IF (initialstep == 0) THEN
!!! if we come from a restart we shall properly compute the average force
!!!      read the average force up to now
         DO i = 1, force_env%meta_env%n_colvar
            fs_section => section_vals_get_subs_vals(force_env%meta_env%metadyn_section, "EXT_LAGRANGE_FS")
            CALL section_vals_get(fs_section, explicit=explicit)
            IF (explicit) THEN
               CALL section_vals_val_get(fs_section, "_DEFAULT_KEYWORD_", &
                                         i_rep_val=i, r_val=rval)
               fz(i) = rval*rand2skip
            END IF
         END DO
 
         CALL hmcsampler(globenv, force_env, mcaverages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, &
                         fz, zbuff, nskip=rand2skip)
         CALL cp_iterate(logger%iter_info, last=.false., iter_nr=0)
         CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip + nmccycles)
         CALL write_restart(md_env=md_env, root_section=force_env%root_section)
      END IF
      IF (output_unit > 0) THEN
         WRITE (output_unit, '(a)') "HMC|==== end initial average forces"
      END IF
!     call set_md_env(md_env, init=.FALSE.)
 
      CALL metadyn_write_colvar(force_env)
 
      DO istep = 1, force_env%meta_env%TAMCSteps
         ! Increase counters
         itimes = itimes + 1
         time = time + force_env%meta_env%dt
         IF (output_unit > 0) THEN
            WRITE (output_unit, '(a)') "HMC|==================================="
            WRITE (output_unit, '(a,1x,i0)') "HMC| on z step ", istep
         END IF
         !needed when electric field fields are applied
         IF (ASSOCIATED(force_env%qs_env)) THEN
            force_env%qs_env%sim_time = time
            force_env%qs_env%sim_step = itimes
            force_env%meta_env%time = force_env%qs_env%sim_time
         END IF
 
         CALL cp_iterate(logger%iter_info, last=(istep == force_env%meta_env%TAMCSteps), iter_nr=itimes)
         ! Open possible Shake output units
         simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, "CONSTRAINT_INFO", &
                                                       extension=".shakeLog", log_filename=.false.)
         simpar%lagrange_multipliers = cp_print_key_unit_nr( &
                                       logger, constraint_section, &
                                       "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.false.)
         simpar%dump_lm = btest(cp_print_key_should_output(logger%iter_info, constraint_section, &
                                                           "LAGRANGE_MULTIPLIERS"), cp_p_file)
 
         ! Velocity Verlet Integrator
 
         moves%hmc%attempts = 0
         moves%hmc%successes = 0
         CALL langevinvec(md_env, globenv, mc_env, moves, gmoves, r, &
                          rng_stream_mc, xieta, an, fz, mcaverages, zbuff)
 
         ! Close Shake output if requested...
         CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
                                           constraint_section, "CONSTRAINT_INFO")
         CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
                                           constraint_section, "LAGRANGE_MULTIPLIERS")
         CALL cp_iterate(logger%iter_info, iter_nr=initialstep)
         CALL metadyn_write_colvar(force_env)
         ! Free Energy calculation
!        CALL free_energy_evaluate(md_env,should_stop,free_energy_section)
 
         ![AME:UB] IF (should_stop) EXIT
 
         ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit
         ! Default:
         ! IF so we don't overwrite the restart or append to the trajectory
         ! because the execution could in principle stop inside the SCF where energy
         ! and forces are not converged.
         ! But:
         ! You can force to print the last step (for example if the method used
         ! to compute energy and forces is not SCF based) activating the print_key
         ! MOTION%MD%PRINT%FORCE_LAST.
         CALL external_control(should_stop, "MD", globenv=globenv)
         IF (should_stop) THEN
            CALL cp_iterate(logger%iter_info, last=.true., iter_nr=itimes)
!           CALL md_output(md_env,md_section,force_env%root_section,should_stop)
            EXIT
         END IF
 
!        IF(simpar%ensemble /= reftraj_ensemble) THEN
!           CALL md_energy(md_env, md_ener)
!           CALL temperature_control(simpar, md_env, md_ener, force_env, logger)
!           CALL comvel_control(md_ener, force_env, md_section, logger)
!           CALL angvel_control(md_ener, force_env, md_section, logger)
!        ELSE
!           CALL md_ener_reftraj(md_env, md_ener)
!        END IF
 
         time_iter_stop = m_walltime()
         used_time = time_iter_stop - time_iter_start
         time_iter_start = time_iter_stop
 
!!!!! this writes the restart...
!         CALL md_output(md_env,md_section,force_env%root_section,should_stop)
 
!        IF(simpar%ensemble == reftraj_ensemble ) THEN
!           CALL write_output_reftraj(md_env)
!        END IF
 
         IF (output_unit > 0) THEN
            WRITE (output_unit, '(a,1x,i0)') "HMC| end z step ", istep
            WRITE (output_unit, '(a)') "HMC|==================================="
         END IF
      END DO
      CALL cp_iterate(logger%iter_info, last=.true., iter_nr=itimes)
      force_env%qs_env%sim_time = 0.0_dp
      force_env%qs_env%sim_step = 0
      rand2skip = rand2skip + nmccycles*force_env%meta_env%TAMCSteps
      IF (initialstep == 0) rand2skip = rand2skip + nmccycles
      CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip)
 
      CALL write_restart(md_env=md_env, root_section=force_env%root_section)
! if we need the final kinetic energy for Hybrid Monte Carlo
!     hmc_ekin%final_ekin=md_ener%ekin
 
      ! Remove the iteration level
      CALL cp_rm_iter_level(logger%iter_info, "MD")
 
      ! Deallocate Thermostats and Barostats
      CALL release_simpar_type(simpar)
 
      CALL md_env_release(md_env)
      DEALLOCATE (md_env)
      ! Clean restartable sections..
      IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
!     IF (para_env%is_source()) THEN
!     ENDIF
      CALL mc_env_release(mc_env)
      DEALLOCATE (mc_env)
      DEALLOCATE (mc_par)
      CALL mc_moves_release(moves)
      CALL mc_moves_release(gmoves)
      DEALLOCATE (r)
!     DEALLOCATE(r_old)
      DEALLOCATE (xieta)
      DEALLOCATE (an)
      DEALLOCATE (fz)
      DEALLOCATE (zbuff)
      CALL mc_averages_release(mcaverages)
      CALL timestop(handle)
 
   END SUBROUTINE qs_tamc
 
! **************************************************************************************************
!> \brief Propagates velocities for z half a step
!> \param force_env ...
!> \param An ...
!> \author Alin M Elena
!> \note   Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
! **************************************************************************************************
   SUBROUTINE tamc_velocities_colvar(force_env, An)
      TYPE(force_env_type), POINTER                      :: force_env
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: an
 
      CHARACTER(len=*), PARAMETER :: routinen = 'tamc_velocities_colvar'
 
      INTEGER                                            :: handle, i_c
      REAL(kind=dp)                                      :: dt, fft, sigma
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(meta_env_type), POINTER                       :: meta_env
      TYPE(metavar_type), POINTER                        :: cv
 
      NULLIFY (logger, meta_env, cv)
      meta_env => force_env%meta_env
      CALL timeset(routinen, handle)
      logger => cp_get_default_logger()
      ! Add citation
      IF (meta_env%langevin) CALL cite_reference(vandencic2006)
      dt = meta_env%dt
 
      ! Evolve Velocities
      meta_env%epot_walls = 0.0_dp
      DO i_c = 1, meta_env%n_colvar
         cv => meta_env%metavar(i_c)
         fft = cv%ff_s + cv%ff_hills
         sigma = sqrt((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
         cv%vvp = cv%vvp + 0.5_dp*dt*(fft/cv%mass - cv%gamma*cv%vvp)*(1.0_dp - 0.25_dp*dt*cv%gamma) + an(i_c)
         meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
      END DO
      CALL timestop(handle)
   END SUBROUTINE tamc_velocities_colvar
 
! **************************************************************************************************
!> \brief propagates z one step
!> \param force_env ...
!> \param xieta ...
!> \author Alin M Elena
!> \note  Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
! **************************************************************************************************
   SUBROUTINE tamc_position_colvar(force_env, xieta)
      TYPE(force_env_type), POINTER                      :: force_env
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: xieta
 
      CHARACTER(len=*), PARAMETER :: routinen = 'tamc_position_colvar'
 
      INTEGER                                            :: handle, i_c
      REAL(kind=dp)                                      :: dt, sigma
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(meta_env_type), POINTER                       :: meta_env
      TYPE(metavar_type), POINTER                        :: cv
 
      NULLIFY (logger, meta_env, cv)
      meta_env => force_env%meta_env
!     IF (.NOT.ASSOCIATED(meta_env)) RETURN
 
      CALL timeset(routinen, handle)
      logger => cp_get_default_logger()
 
      ! Add citation
      IF (meta_env%langevin) CALL cite_reference(vandencic2006)
      dt = meta_env%dt
 
      ! Update of ss0
      DO i_c = 1, meta_env%n_colvar
         cv => meta_env%metavar(i_c)
         sigma = sqrt((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
!        cv%ss0 =cv%ss0 +dt*cv%vvp
         cv%ss0 = cv%ss0 + dt*cv%vvp + dt*sqrt(dt/12.0_dp)*sigma*xieta(i_c + meta_env%n_colvar)
         IF (cv%periodic) THEN
            ! A periodic COLVAR is always within [-pi,pi]
            cv%ss0 = sign(1.0_dp, asin(sin(cv%ss0)))*acos(cos(cv%ss0))
         END IF
      END DO
      CALL timestop(handle)
 
   END SUBROUTINE tamc_position_colvar
 
! **************************************************************************************************
!> \brief Computes forces on z
!> #details also can be used to get the potenzial evergy of z
!> \param force_env ...
!> \param zpot ...
!> \author Alin M Elena
! **************************************************************************************************
   SUBROUTINE tamc_force(force_env, zpot)
      TYPE(force_env_type), POINTER                      :: force_env
      REAL(kind=dp), INTENT(inout), OPTIONAL             :: zpot
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'tamc_force'
 
      INTEGER                                            :: handle, i, i_c, icolvar, ii
      LOGICAL                                            :: explicit
      REAL(kind=dp)                                      :: diff_ss, dt, rval
      TYPE(colvar_p_type), DIMENSION(:), POINTER         :: colvar_p
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(meta_env_type), POINTER                       :: meta_env
      TYPE(metavar_type), POINTER                        :: cv
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(section_vals_type), POINTER                   :: ss0_section, ss_section, vvp_section
 
      NULLIFY (logger, meta_env)
      meta_env => force_env%meta_env
!     IF (.NOT.ASSOCIATED(meta_env)) RETURN
 
      CALL timeset(routinen, handle)
      logger => cp_get_default_logger()
      NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section, ss_section)
      CALL force_env_get(force_env, subsys=subsys)
 
      dt = meta_env%dt
      IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1
      ! compute ss and the derivative of ss with respect to the atomic positions
      DO i_c = 1, meta_env%n_colvar
         cv => meta_env%metavar(i_c)
         icolvar = cv%icolvar
         CALL colvar_eval_glob_f(icolvar, force_env)
         cv%ss = subsys%colvar_p(icolvar)%colvar%ss
         ! Restart for Extended Lagrangian Metadynamics
         IF (meta_env%restart) THEN
            ! Initialize the position of the collective variable in the extended lagrange
            ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
            CALL section_vals_get(ss0_section, explicit=explicit)
            IF (explicit) THEN
               CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", &
                                         i_rep_val=i_c, r_val=rval)
               cv%ss0 = rval
            ELSE
               cv%ss0 = cv%ss
            END IF
            vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
            CALL section_vals_get(vvp_section, explicit=explicit)
            IF (explicit) THEN
               CALL setup_velocities_z(force_env)
               CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", &
                                         i_rep_val=i_c, r_val=rval)
               cv%vvp = rval
            ELSE
               CALL setup_velocities_z(force_env)
            END IF
            ss_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS")
            CALL section_vals_get(ss_section, explicit=explicit)
            IF (explicit) THEN
               CALL section_vals_val_get(ss_section, "_DEFAULT_KEYWORD_", &
                                         i_rep_val=i_c, r_val=rval)
               cv%ss = rval
            END IF
         END IF
         !
      END DO
      ! forces on the atoms
      NULLIFY (particles)
      CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
                         particles=particles)
 
      meta_env%restart = .false.
      meta_env%epot_s = 0.0_dp
      meta_env%epot_walls = 0.0_dp
      DO i_c = 1, meta_env%n_colvar
         cv => meta_env%metavar(i_c)
         diff_ss = cv%ss - cv%ss0
         IF (cv%periodic) THEN
            ! The difference of a periodic COLVAR is always within [-pi,pi]
            diff_ss = sign(1.0_dp, asin(sin(diff_ss)))*acos(cos(diff_ss))
         END IF
         cv%epot_s = 0.5_dp*cv%lambda*diff_ss*diff_ss
         cv%ff_s = cv%lambda*(diff_ss)
         meta_env%epot_s = meta_env%epot_s + cv%epot_s
         icolvar = cv%icolvar
 
         DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
            i = colvar_p(icolvar)%colvar%i_atom(ii)
            particles%els(i)%f = particles%els(i)%f - cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:, ii)
         END DO
 
      END DO
      IF (PRESENT(zpot)) zpot = meta_env%epot_s
      CALL fix_atom_control(force_env)
 
      CALL timestop(handle)
   END SUBROUTINE tamc_force
 
! **************************************************************************************************
!> \brief propagates one time step both z systems and samples the x system
!> \param md_env ...
!> \param globenv ...
!> \param mc_env ...
!> \param moves ...
!> \param gmoves ...
!> \param r ...
!> \param rng_stream_mc ...
!> \param xieta ...
!> \param An ...
!> \param fz ...
!> \param averages ...
!> \param zbuff ...
!> \author Alin M Elena
! **************************************************************************************************
   SUBROUTINE langevinvec(md_env, globenv, mc_env, moves, gmoves, r, &
                          rng_stream_mc, xieta, An, fz, averages, zbuff)
 
      TYPE(md_environment_type), POINTER                 :: md_env
      TYPE(global_environment_type), POINTER             :: globenv
      TYPE(mc_environment_type), POINTER                 :: mc_env
      TYPE(mc_moves_type), POINTER                       :: moves, gmoves
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: r
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream_mc
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: xieta, an, fz
      TYPE(mc_averages_type), INTENT(INOUT), POINTER     :: averages
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: zbuff
 
      INTEGER                                            :: iprint, ivar, nparticle, nparticle_kind, &
                                                            nstep, output_unit
      REAL(kind=dp)                                      :: dt, gamma, mass, sigma
      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(global_constraint_type), POINTER              :: gci
      TYPE(mc_simpar_type), POINTER                      :: mc_par
      TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
      TYPE(molecule_list_type), POINTER                  :: molecules
      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(simpar_type), POINTER                         :: simpar
      TYPE(virial_type), POINTER                         :: virial
 
      NULLIFY (logger, mc_par)
      logger => cp_get_default_logger()
      output_unit = cp_logger_get_default_io_unit(logger)
 
! quantitites to be nullified for the get_md_env
      NULLIFY (simpar, force_env, para_env)
! quantities to be nullified for the force_env_get environment
      NULLIFY (subsys, cell)
! quantitites to be nullified for the cp_subsys_get
      NULLIFY (atomic_kinds, local_particles, particles, local_molecules, molecules, molecule_kinds, gci)
 
      CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
                      para_env=para_env)
      CALL get_mc_env(mc_env, mc_par=mc_par)
      CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint)
 
      dt = simpar%dt
      CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
 
!!!! this bit should vanish once I understand what the hell is with it
 
!     ! Do some checks on coordinates and box
      CALL apply_qmmm_walls_reflective(force_env)
 
      CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
                         particles=particles, local_molecules=local_molecules, molecules=molecules, &
                         molecule_kinds=molecule_kinds, gci=gci, virial=virial)
 
      nparticle_kind = atomic_kinds%n_els
      atomic_kind_set => atomic_kinds%els
      molecule_kind_set => molecule_kinds%els
 
      nparticle = particles%n_els
      particle_set => particles%els
      molecule_set => molecules%els
      cpassert(ASSOCIATED(force_env%meta_env))
      cpassert(force_env%meta_env%langevin)
      !    *** Velocity Verlet for Langevin *** v(t)--> v(t+1/2)
      !!!!!! noise xi is in the first half, eta in the second half
      DO ivar = 1, force_env%meta_env%n_colvar
         xieta(ivar) = force_env%meta_env%rng(ivar)%next()
         xieta(ivar + force_env%meta_env%n_colvar) = force_env%meta_env%rng(ivar)%next()
         gamma = force_env%meta_env%metavar(ivar)%gamma
         mass = force_env%meta_env%metavar(ivar)%mass
         sigma = sqrt((force_env%meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*gamma/mass)
         an(ivar) = 0.5_dp*sqrt(dt)*sigma*(xieta(ivar)*(1.0_dp - 0.5_dp*dt*gamma) - &
                                           dt*gamma*xieta(ivar + force_env%meta_env%n_colvar)/sqrt(12.0_dp))
      END DO
!    *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
      CALL tamc_velocities_colvar(force_env, an)
!    *** Velocity Verlet for Langevin S(t)->S(t+1)
      CALL tamc_position_colvar(force_env, xieta)
!!!!! start zHMC sampler
      CALL hmcsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, fz, zbuff)
 
!     CALL final_mc_write(mc_par,tmp_moves,&
!                output_unit,energy_check,&
!                initial_energy,final_energy,&
!                averages)
 
!!!!!!!!!!!!!!!!!!!! end zHMC sampler
      !    *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t+1)
      CALL tamc_velocities_colvar(force_env, an)
!       CALL virial_evaluate ( atomic_kind_set, particle_set,  &
!          local_particles, virial, para_env)
 
   END SUBROUTINE langevinvec
 
! **************************************************************************************************
!> \brief Driver routin for the canonical sampler using modified HMC
!> \param globenv ...
!> \param force_env ...
!> \param averages ...
!> \param r ...
!> \param mc_par ...
!> \param moves ...
!> \param gmoves ...
!> \param rng_stream_mc ...
!> \param output_unit ...
!> \param fz ...
!> \param zbuff ...
!> \param nskip ...
!> \author Alin M Elena
!> \note at the end of this routine %ff_s will contain mean force
! **************************************************************************************************
 
   SUBROUTINE hmcsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, &
                         fz, zbuff, nskip)
      TYPE(global_environment_type), POINTER             :: globenv
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(mc_averages_type), POINTER                    :: averages
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: r
      TYPE(mc_simpar_type), POINTER                      :: mc_par
      TYPE(mc_moves_type), POINTER                       :: moves, gmoves
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream_mc
      INTEGER, INTENT(IN)                                :: output_unit
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: fz, zbuff
      INTEGER, INTENT(IN), OPTIONAL                      :: nskip
 
      INTEGER                                            :: i, iprint, ishift, it1, j, nsamples, &
                                                            nstep
      REAL(kind=dp)                                      :: energy_check, old_epx, old_epz, t1
      TYPE(meta_env_type), POINTER                       :: meta_env_saved
 
      IF (PRESENT(nskip)) THEN
         nsamples = nskip
         ishift = nskip
      ELSE
         nsamples = 0
         fz = 0.0_dp
         ishift = 0
      END IF
!      lrestart = .false.
!      if (present(logger) .and. present(iter)) THEN
!       lrestart=.true.
!      ENDIF
      CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint)
      meta_env_saved => force_env%meta_env
      NULLIFY (force_env%meta_env)
      CALL force_env_get(force_env, potential_energy=old_epx)
      force_env%meta_env => meta_env_saved
 
      old_epz = force_env%meta_env%epot_s
!!! average energy will be wrong on restarts
      averages%ave_energy = 0.0_dp
      t1 = force_env%qs_env%sim_time
      it1 = force_env%qs_env%sim_step
      IF (output_unit > 0) THEN
         WRITE (output_unit, '(a,l4)') "HMC|restart? ", force_env%meta_env%restart
         WRITE (output_unit, '(a,3(f16.8,1x))') &
            "HMC|Ep, Epx, Epz ", old_epx + force_env%meta_env%epot_s, old_epx, force_env%meta_env%epot_s
         WRITE (output_unit, '(a)') "#HMC| No | z.. | theta.. | ff_z... | ff_z/n |"
      END IF
      DO i = 1, nstep
         IF (mod(i, iprint) == 0 .AND. (output_unit > 0)) THEN
            WRITE (output_unit, '(a,1x,i0)') "HMC|========== On Monte Carlo cycle ", i + ishift
            WRITE (output_unit, '(a)') "HMC| Attempting a minitrajectory move"
            WRITE (output_unit, '(a,1x,i0)') "HMC| start mini-trajectory", i + ishift
            WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|0 ", i + ishift
            DO j = 1, force_env%meta_env%n_colvar
               WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, &
                  force_env%meta_env%metavar(j)%ss, &
                  force_env%meta_env%metavar(j)%ff_s !,fz(j)/real(i+ishift,dp)
            END DO
            WRITE (output_unit, *)
         END IF
 
         CALL mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, energy_check, &
                          r, output_unit, rng_stream_mc, zbuff)
         ! check averages...
         ! force average for z needed too...
         averages%ave_energy = averages%ave_energy*real(i - 1, dp)/real(i, dp) + &
                               old_epx/real(i, dp)
         DO j = 1, force_env%meta_env%n_colvar
            fz(j) = fz(j) + force_env%meta_env%metavar(j)%ff_s
         END DO
         IF (output_unit > 0) THEN
            WRITE (output_unit, '(a,1x,i0)') "HMC|end mini-trajectory", i + ishift
!!!!!!!! this prints z and theta(x) --ss0,ss-- needed to determine an acceptable k then
            !  the instanteneous force and some instanteneous average for force
            WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|1 ", i + ishift
            DO j = 1, force_env%meta_env%n_colvar
               WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, &
                  force_env%meta_env%metavar(j)%ss, &
                  force_env%meta_env%metavar(j)%ff_s, fz(j)/real(i + ishift, dp)
            END DO
            WRITE (output_unit, *)
         END IF
         nsamples = nsamples + 1
         IF (mod(i, iprint) == 0 .AND. (output_unit > 0)) THEN
            WRITE (output_unit, '(a,f16.8)') "HMC| Running average for potential energy ", averages%ave_energy
            WRITE (output_unit, '(a,1x,i0)') "HMC|======== End Monte Carlo cycle ", i + ishift
         END IF
!         IF (lrestart) THEN
!           k=nstep/5
!           IF(MOD(i,k) == 0) THEN
!              force_env%qs_env%sim_time=t1
!              force_env%qs_env%sim_step=it1
!              DO j=1,force_env%meta_env%n_colvar
!                force_env%meta_env%metavar(j)%ff_s=fz(j)/real(i+ishift,dp)
!              ENDDO
! !              CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=-1)
!              CALL section_vals_val_set(mcsec,"RANDOMTOSKIP",i_val=i+ishift)
!              CALL write_restart(md_env=mdenv,root_section=force_env%root_section)
! !              CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=iter)
!           ENDIF
!         ENDIF
      END DO
      force_env%qs_env%sim_time = t1
      force_env%qs_env%sim_step = it1
      IF (output_unit > 0) THEN
         WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| local acceptance ratio: ", moves%hmc%successes, "/", &
            moves%hmc%attempts, "=", real(moves%hmc%successes, dp)/real(moves%hmc%attempts, dp)
         WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| global acceptance ratio: ", gmoves%hmc%successes, "/", &
            gmoves%hmc%attempts, "=", real(gmoves%hmc%successes, dp)/real(gmoves%hmc%attempts, dp)
      END IF
      !average force
      DO j = 1, force_env%meta_env%n_colvar
         force_env%meta_env%metavar(j)%ff_s = fz(j)/nsamples
      END DO
   END SUBROUTINE hmcsampler
 
! **************************************************************************************************
!> \brief performs a hybrid Monte Carlo move
!> \param mc_par ...
!> \param force_env ...
!> \param globenv ...
!> \param moves ...
!> \param gmoves ...
!> \param old_epx ...
!> \param old_epz ...
!> \param energy_check ...
!> \param r ...
!> \param output_unit ...
!> \param rng_stream ...
!> \param zbuff ...
!> \author Alin M Elena
!> \note It runs a NVE trajectory, followed by localisation and accepts rejects
!> using the biased Hamiltonian, rather than the traditional guiding Hamiltonian
! **************************************************************************************************
   SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, &
                          energy_check, r, output_unit, rng_stream, zbuff)
 
      TYPE(mc_simpar_type), POINTER                      :: mc_par
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(global_environment_type), POINTER             :: globenv
      TYPE(mc_moves_type), POINTER                       :: moves, gmoves
      REAL(kind=dp), INTENT(INOUT)                       :: old_epx, old_epz, energy_check
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: r
      INTEGER, INTENT(IN)                                :: output_unit
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: zbuff
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'mc_hmc_move'
 
      INTEGER                                            :: handle, iatom, j, natoms, source
      LOGICAL                                            :: ionode, localise
      REAL(kind=dp)                                      :: beta, energy_term, exp_max_val, &
                                                            exp_min_val, new_energy, new_epx, &
                                                            new_epz, rand, value, w
      TYPE(cp_subsys_type), POINTER                      :: oldsys
      TYPE(mc_ekin_type), POINTER                        :: hmc_ekin
      TYPE(meta_env_type), POINTER                       :: meta_env_saved
      TYPE(mp_comm_type)                                 :: group
      TYPE(particle_list_type), POINTER                  :: particles_set
      TYPE(section_vals_type), POINTER                   :: dft_section, input
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(force_env%qs_env, input=input)
      dft_section => section_vals_get_subs_vals(input, "DFT")
 
! get a bunch of stuff from mc_par
      CALL get_mc_par(mc_par, ionode=ionode, &
                      beta=beta, exp_max_val=exp_max_val, &
                      exp_min_val=exp_min_val, source=source, group=group)
 
! nullify some pointers
!       NULLIFY(particles_set,oldsys,hmc_ekin)
      NULLIFY (particles_set, oldsys, meta_env_saved, hmc_ekin)
      ! now let's grab the particle positions
      CALL force_env_get(force_env, subsys=oldsys)
      CALL cp_subsys_get(oldsys, particles=particles_set)
      natoms = SIZE(particles_set%els)
! do some allocation
 
      ALLOCATE (hmc_ekin)
 
! record the attempt
      moves%hmc%attempts = moves%hmc%attempts + 1
      gmoves%hmc%attempts = gmoves%hmc%attempts + 1
 
! save the old coordinates just in case we need to go back
      DO iatom = 1, natoms
         r(1:3, iatom) = particles_set%els(iatom)%r(1:3)
      END DO
      localise = .true.
! the same for collective variables data should be made,ss first half and ff_s the last half
      DO j = 1, force_env%meta_env%n_colvar
         zbuff(j) = force_env%meta_env%metavar(j)%ss
         zbuff(j + force_env%meta_env%n_colvar) = force_env%meta_env%metavar(j)%ff_s
         IF ((oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == hbp_colvar_id) .OR. &
             (oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == wc_colvar_id)) THEN
            localise = .false.
         END IF
      END DO
 
! now run the MD simulation
      meta_env_saved => force_env%meta_env
      NULLIFY (force_env%meta_env)
      force_env%qs_env%sim_time = 0.0_dp
      force_env%qs_env%sim_step = 0
      IF (.NOT. localise) THEN
         CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.false.)
      END IF
      CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin)
      IF (.NOT. localise) THEN
         CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.true.)
         CALL scf_post_calculation_gpw(force_env%qs_env)
      END IF
 
      CALL force_env_get(force_env, potential_energy=new_epx)
 
      force_env%meta_env => meta_env_saved
      CALL tamc_force(force_env, zpot=new_epz)
      new_energy = new_epx + new_epz
      IF (output_unit > 0) THEN
         WRITE (output_unit, '(a,4(f16.8,1x))') &
            "HMC|old Ep, Ekx, Epz, Epx ", old_epx + old_epz, hmc_ekin%initial_ekin, old_epz, old_epx
         WRITE (output_unit, '(a,4(f16.8,1x))') "HMC|new Ep, Ekx, Epz, Epx ", new_energy, hmc_ekin%final_ekin, new_epz, new_epx
      END IF
      energy_term = new_energy - old_epx - old_epz + hmc_ekin%final_ekin - hmc_ekin%initial_ekin
 
      value = -beta*(energy_term)
! to prevent overflows
      IF (value > exp_max_val) THEN
         w = 10.0_dp
      ELSEIF (value < exp_min_val) THEN
         w = 0.0_dp
      ELSE
         w = exp(value)
      END IF
 
      rand = rng_stream%next()
      IF (rand < w) THEN
! accept the move
         moves%hmc%successes = moves%hmc%successes + 1
         gmoves%hmc%successes = gmoves%hmc%successes + 1
! update energies
         energy_check = energy_check + (new_energy - old_epx - old_epz)
         old_epx = new_epx
         old_epz = new_epz
      ELSE
! reset the cell and particle positions
         DO iatom = 1, natoms
            particles_set%els(iatom)%r(1:3) = r(1:3, iatom)
         END DO
         DO j = 1, force_env%meta_env%n_colvar
            force_env%meta_env%metavar(j)%ss = zbuff(j)
            force_env%meta_env%metavar(j)%ff_s = zbuff(j + force_env%meta_env%n_colvar)
         END DO
 
      END IF
 
      DEALLOCATE (hmc_ekin)
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE mc_hmc_move
 
! **************************************************************************************************
!> \brief ...
!> \param force_env ...
! **************************************************************************************************
   SUBROUTINE metadyn_write_colvar_header(force_env)
      TYPE(force_env_type), POINTER                      :: force_env
 
      CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_write_colvar_header'
 
      CHARACTER(len=100)                                 :: aux, fmt
      CHARACTER(len=255)                                 :: label1, label2, label3, label4, label5, &
                                                            label6
      INTEGER                                            :: handle, i, iw, m
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(meta_env_type), POINTER                       :: meta_env
 
      NULLIFY (logger, meta_env)
      meta_env => force_env%meta_env
      IF (.NOT. ASSOCIATED(meta_env)) RETURN
 
      CALL timeset(routinen, handle)
      logger => cp_get_default_logger()
 
      iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
                                "PRINT%COLVAR", extension=".metadynLog")
      IF (iw > 0) THEN
         label1 = ""
         label2 = ""
         label3 = ""
         label4 = ""
         label5 = ""
         label6 = ""
         DO i = 1, meta_env%n_colvar
            WRITE (aux, '(a,i0)') "z_", i
            label1 = trim(label1)//trim(aux)
            m = 15*i - len_trim(label1) - 1
            label1 = trim(label1)//repeat(" ", m)//"|"
            WRITE (aux, '(a,i0)') "Theta_", i
            label2 = trim(label2)//trim(aux)
            m = 15*i - len_trim(label2) - 1
            label2 = trim(label2)//repeat(" ", m)//"|"
            WRITE (aux, '(a,i0)') "F_z", i
            label3 = trim(label3)//trim(aux)
            m = 15*i - len_trim(label3) - 1
            label3 = trim(label3)//repeat(" ", m)//"|"
            WRITE (aux, '(a,i0)') "F_h", i
            label4 = trim(label4)//trim(aux)
            m = 15*i - len_trim(label4) - 1
            label4 = trim(label4)//repeat(" ", m)//"|"
            WRITE (aux, '(a,i0)') "F_w", i
            label5 = trim(label5)//trim(aux)
            m = 15*i - len_trim(label5) - 1
            label5 = trim(label5)//repeat(" ", m)//"|"
            WRITE (aux, '(a,i0)') "v_z", i
            label6 = trim(label6)//trim(aux)
            m = 15*i - len_trim(label6) - 1
            label6 = trim(label6)//repeat(" ", m)//"|"
         END DO
         WRITE (fmt, '("(a17,6a",i0 ,",4a15)")') meta_env%n_colvar*15
         WRITE (iw, trim(fmt)) "#Time[fs] |", &
            trim(label1), &
            trim(label2), &
            trim(label3), &
            trim(label4), &
            trim(label5), &
            trim(label6), &
            "Epot_z |", &
            "Ene hills |", &
            "Epot walls |", &
            "Temperature |"
 
      END IF
      CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
                                        "PRINT%COLVAR")
 
      CALL timestop(handle)
 
   END SUBROUTINE metadyn_write_colvar_header
 
! **************************************************************************************************
!> \brief ...
!> \param force_env ...
! **************************************************************************************************
   SUBROUTINE metadyn_write_colvar(force_env)
      TYPE(force_env_type), POINTER                      :: force_env
 
      CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_write_colvar'
 
      INTEGER                                            :: handle, i, i_c, iw
      REAL(kind=dp)                                      :: temp
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(meta_env_type), POINTER                       :: meta_env
      TYPE(metavar_type), POINTER                        :: cv
 
      NULLIFY (logger, meta_env, cv)
      meta_env => force_env%meta_env
      IF (.NOT. ASSOCIATED(meta_env)) RETURN
 
      CALL timeset(routinen, handle)
      logger => cp_get_default_logger()
 
      IF (meta_env%langevin) THEN
         meta_env%ekin_s = 0.0_dp
!        meta_env%epot_s = 0.0_dp
         DO i_c = 1, meta_env%n_colvar
            cv => meta_env%metavar(i_c)
            meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
         END DO
      END IF
 
      ! write COLVAR file
      iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
                                "PRINT%COLVAR", extension=".metadynLog")
      IF (iw > 0) THEN
         IF (meta_env%extended_lagrange) THEN
            WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, &
               (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), &
               meta_env%epot_s, &
               meta_env%hills_env%energy, &
               meta_env%epot_walls, &
               (meta_env%ekin_s)*2.0_dp/(real(meta_env%n_colvar, kind=dp))*kelvin
         ELSE
            WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, &
               (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
               (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
               meta_env%hills_env%energy, &
               meta_env%epot_walls
         END IF
      END IF
      CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
                                        "PRINT%COLVAR")
      ! Temperature for COLVAR
      IF (meta_env%extended_lagrange) THEN
         temp = meta_env%ekin_s*2.0_dp/(real(meta_env%n_colvar, kind=dp))*kelvin
         meta_env%avg_temp = (meta_env%avg_temp*real(meta_env%n_steps, kind=dp) + &
                              temp)/real(meta_env%n_steps + 1, kind=dp)
         iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
                                   "PRINT%TEMPERATURE_COLVAR", extension=".metadynLog")
         IF (iw > 0) THEN
            WRITE (iw, '(T2,79("-"))')
            WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', &
               temp, meta_env%avg_temp
            WRITE (iw, '(T2,79("-"))')
         END IF
         CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
                                           "PRINT%TEMPERATURE_COLVAR")
      END IF
      CALL timestop(handle)
 
   END SUBROUTINE metadyn_write_colvar
 
! **************************************************************************************************
!> \brief ...
!> \param force_env ...
! **************************************************************************************************
   SUBROUTINE setup_velocities_z(force_env)
      TYPE(force_env_type), POINTER                      :: force_env
 
      INTEGER                                            :: i_c
      REAL(kind=dp)                                      :: ekin_w, fac_t
      TYPE(meta_env_type), POINTER                       :: meta_env
      TYPE(metavar_type), POINTER                        :: cv
 
      NULLIFY (meta_env)
      meta_env => force_env%meta_env
      meta_env%ekin_s = 0.0_dp
      DO i_c = 1, meta_env%n_colvar
         cv => meta_env%metavar(i_c)
         cv%vvp = force_env%globenv%gaussian_rng_stream%next()
         meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
      END DO
      ekin_w = 0.5_dp*meta_env%temp_wanted*real(meta_env%n_colvar, kind=dp)
      fac_t = sqrt(ekin_w/max(meta_env%ekin_s, 1.0e-8_dp))
      DO i_c = 1, meta_env%n_colvar
         cv => meta_env%metavar(i_c)
         cv%vvp = cv%vvp*fac_t
      END DO
   END SUBROUTINE setup_velocities_z
END MODULE tamc_run