LCOV - code coverage report
Current view: top level - src/motion - pint_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 1052 1226 85.8 %
Date: 2023-03-30 11:55:16 Functions: 19 21 90.5 %

          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  Methods to performs a path integral run
      10             : !> \author fawzi
      11             : !> \par History
      12             : !>      02.2005 created [fawzi]
      13             : !>           11.2006 modified so it might actually work [hforbert]
      14             : !>           10.2015 added RPMD propagator
      15             : !>           10.2015 added exact harmonic integrator [Felix Uhl]
      16             : !> \note   quick & dirty rewrite of my python program
      17             : ! **************************************************************************************************
      18             : MODULE pint_methods
      19             : 
      20             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      21             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      22             :                                               get_atomic_kind
      23             :    USE bibliography,                    ONLY: Brieuc2016,&
      24             :                                               Ceriotti2010,&
      25             :                                               Ceriotti2012,&
      26             :                                               cite_reference
      27             :    USE cell_types,                      ONLY: cell_type
      28             :    USE constraint,                      ONLY: rattle_control,&
      29             :                                               shake_control,&
      30             :                                               shake_update_targets
      31             :    USE constraint_util,                 ONLY: getold
      32             :    USE cp_external_control,             ONLY: external_control
      33             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      34             :                                               cp_logger_get_default_io_unit,&
      35             :                                               cp_logger_type,&
      36             :                                               cp_to_string
      37             :    USE cp_output_handling,              ONLY: cp_add_iter_level,&
      38             :                                               cp_iterate,&
      39             :                                               cp_p_file,&
      40             :                                               cp_print_key_finished_output,&
      41             :                                               cp_print_key_should_output,&
      42             :                                               cp_print_key_unit_nr,&
      43             :                                               cp_rm_iter_level
      44             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      45             :                                               cp_subsys_type
      46             :    USE cp_units,                        ONLY: cp_unit_from_cp2k,&
      47             :                                               cp_unit_to_cp2k
      48             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      49             :    USE f77_interface,                   ONLY: f_env_add_defaults,&
      50             :                                               f_env_rm_defaults,&
      51             :                                               f_env_type
      52             :    USE force_env_types,                 ONLY: force_env_get
      53             :    USE gle_system_dynamics,             ONLY: gle_cholesky_stab,&
      54             :                                               gle_matrix_exp,&
      55             :                                               restart_gle
      56             :    USE gle_system_types,                ONLY: gle_dealloc,&
      57             :                                               gle_init,&
      58             :                                               gle_thermo_create
      59             :    USE global_types,                    ONLY: global_environment_type
      60             :    USE helium_interactions,             ONLY: helium_intpot_scan
      61             :    USE helium_io,                       ONLY: helium_write_cubefile
      62             :    USE helium_methods,                  ONLY: helium_create,&
      63             :                                               helium_init,&
      64             :                                               helium_release
      65             :    USE helium_sampling,                 ONLY: helium_do_run,&
      66             :                                               helium_step
      67             :    USE helium_types,                    ONLY: helium_solvent_p_type
      68             :    USE input_constants,                 ONLY: integrate_exact,&
      69             :                                               integrate_numeric,&
      70             :                                               propagator_cmd,&
      71             :                                               propagator_rpmd,&
      72             :                                               transformation_normal,&
      73             :                                               transformation_stage
      74             :    USE input_cp2k_restarts,             ONLY: write_restart
      75             :    USE input_section_types,             ONLY: &
      76             :         section_type, section_vals_get, section_vals_get_subs_vals, section_vals_release, &
      77             :         section_vals_retain, section_vals_type, section_vals_val_get, section_vals_val_set, &
      78             :         section_vals_val_unset
      79             :    USE kinds,                           ONLY: default_path_length,&
      80             :                                               default_string_length,&
      81             :                                               dp
      82             :    USE machine,                         ONLY: m_flush,&
      83             :                                               m_walltime
      84             :    USE mathconstants,                   ONLY: twopi
      85             :    USE mathlib,                         ONLY: gcd
      86             :    USE message_passing,                 ONLY: mp_comm_self,&
      87             :                                               mp_para_env_type
      88             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      89             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      90             :    USE molecule_list_types,             ONLY: molecule_list_type
      91             :    USE molecule_types,                  ONLY: global_constraint_type,&
      92             :                                               molecule_type
      93             :    USE parallel_rng_types,              ONLY: GAUSSIAN,&
      94             :                                               rng_stream_type
      95             :    USE particle_list_types,             ONLY: particle_list_type
      96             :    USE particle_types,                  ONLY: particle_type
      97             :    USE pint_gle,                        ONLY: pint_calc_gle_energy,&
      98             :                                               pint_gle_init,&
      99             :                                               pint_gle_step
     100             :    USE pint_io,                         ONLY: pint_write_action,&
     101             :                                               pint_write_centroids,&
     102             :                                               pint_write_com,&
     103             :                                               pint_write_ener,&
     104             :                                               pint_write_line,&
     105             :                                               pint_write_rgyr,&
     106             :                                               pint_write_step_info,&
     107             :                                               pint_write_trajectory
     108             :    USE pint_normalmode,                 ONLY: normalmode_calc_uf_h,&
     109             :                                               normalmode_env_create,&
     110             :                                               normalmode_init_masses,&
     111             :                                               normalmode_release
     112             :    USE pint_piglet,                     ONLY: pint_calc_piglet_energy,&
     113             :                                               pint_piglet_create,&
     114             :                                               pint_piglet_init,&
     115             :                                               pint_piglet_release,&
     116             :                                               pint_piglet_step
     117             :    USE pint_pile,                       ONLY: pint_calc_pile_energy,&
     118             :                                               pint_pile_init,&
     119             :                                               pint_pile_release,&
     120             :                                               pint_pile_step
     121             :    USE pint_public,                     ONLY: pint_levy_walk
     122             :    USE pint_qtb,                        ONLY: pint_calc_qtb_energy,&
     123             :                                               pint_qtb_init,&
     124             :                                               pint_qtb_release,&
     125             :                                               pint_qtb_step
     126             :    USE pint_staging,                    ONLY: staging_calc_uf_h,&
     127             :                                               staging_env_create,&
     128             :                                               staging_init_masses,&
     129             :                                               staging_release
     130             :    USE pint_transformations,            ONLY: pint_f2uf,&
     131             :                                               pint_u2x,&
     132             :                                               pint_x2u
     133             :    USE pint_types,                      ONLY: &
     134             :         e_conserved_id, e_kin_thermo_id, e_kin_virial_id, e_potential_id, pint_env_type, &
     135             :         thermostat_gle, thermostat_none, thermostat_nose, thermostat_piglet, thermostat_pile, &
     136             :         thermostat_qtb
     137             :    USE replica_methods,                 ONLY: rep_env_calc_e_f,&
     138             :                                               rep_env_create
     139             :    USE replica_types,                   ONLY: rep_env_release,&
     140             :                                               replica_env_type
     141             :    USE simpar_types,                    ONLY: create_simpar_type,&
     142             :                                               release_simpar_type
     143             : #include "../base/base_uses.f90"
     144             : 
     145             :    IMPLICIT NONE
     146             :    PRIVATE
     147             : 
     148             :    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
     149             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_methods'
     150             : 
     151             :    PUBLIC :: do_pint_run
     152             : 
     153             : CONTAINS
     154             : 
     155             : ! ***************************************************************************
     156             : !> \brief  Create a path integral environment
     157             : !> \param pint_env ...
     158             : !> \param input ...
     159             : !> \param input_declaration ...
     160             : !> \param para_env ...
     161             : !> \par    History
     162             : !>           Fixed some bugs [hforbert]
     163             : !>           Added normal mode transformation [hforbert]
     164             : !>           10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl]
     165             : !>           10.2018 Added centroid constraints [cschran+rperez]
     166             : !>           10.2021 Added beadwise constraints [lduran]
     167             : !> \author fawzi
     168             : !> \note   Might return an unassociated pointer in parallel on the processors
     169             : !>         that are not needed.
     170             : ! **************************************************************************************************
     171        1350 :    SUBROUTINE pint_create(pint_env, input, input_declaration, para_env)
     172             : 
     173             :       TYPE(pint_env_type), INTENT(OUT)                   :: pint_env
     174             :       TYPE(section_vals_type), POINTER                   :: input
     175             :       TYPE(section_type), POINTER                        :: input_declaration
     176             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     177             : 
     178             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pint_create'
     179             : 
     180             :       CHARACTER(len=2*default_string_length)             :: msg
     181             :       CHARACTER(len=default_path_length)                 :: output_file_name, project_name
     182             :       INTEGER                                            :: handle, iat, ibead, icont, idim, idir, &
     183             :                                                             ierr, ig, itmp, nrep, prep, stat
     184             :       LOGICAL                                            :: explicit, ltmp
     185             :       REAL(kind=dp)                                      :: dt, mass, omega
     186             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     187             :       TYPE(f_env_type), POINTER                          :: f_env
     188             :       TYPE(global_constraint_type), POINTER              :: gci
     189             :       TYPE(particle_list_type), POINTER                  :: particles
     190             :       TYPE(replica_env_type), POINTER                    :: rep_env
     191             :       TYPE(section_vals_type), POINTER :: constraint_section, gle_section, nose_section, &
     192             :          piglet_section, pile_section, pint_section, qtb_section, transform_section
     193             : 
     194          54 :       CALL timeset(routineN, handle)
     195             : 
     196          54 :       NULLIFY (f_env, subsys, particles, nose_section, gle_section, gci)
     197             : 
     198          54 :       CPASSERT(ASSOCIATED(input))
     199          54 :       CPASSERT(input%ref_count > 0)
     200          54 :       NULLIFY (rep_env)
     201          54 :       pint_section => section_vals_get_subs_vals(input, "MOTION%PINT")
     202          54 :       CALL section_vals_val_get(pint_section, "p", i_val=nrep)
     203             :       CALL section_vals_val_get(pint_section, "proc_per_replica", &
     204          54 :                                 i_val=prep)
     205             :       ! Maybe let the user have his/her way as long as prep is
     206             :       ! within the bounds of number of CPUs??
     207          54 :       IF ((prep < 1) .OR. (prep > para_env%num_pe) .OR. &
     208             :           (MOD(prep*nrep, para_env%num_pe) /= 0)) THEN
     209           2 :          prep = para_env%num_pe/gcd(para_env%num_pe, nrep)
     210           2 :          IF (para_env%is_source()) THEN
     211           1 :             WRITE (UNIT=msg, FMT=*) "PINT WARNING: Adjusting number of processors per replica to ", prep
     212          53 :             CPWARN(msg)
     213             :          END IF
     214             :       END IF
     215             : 
     216             :       ! replica_env modifies the global input structure - which is wrong - one
     217             :       ! of the side effects is the inifite adding of the -r-N string to the
     218             :       ! project name and the output file name, which corrupts restart files.
     219             :       ! For now: save the project name and output file name and restore them
     220             :       ! after the rep_env_create has executed - the initialization of the
     221             :       ! replicas will run correctly anyways.
     222             :       ! TODO: modify rep_env so that it behaves better
     223          54 :       CALL section_vals_val_get(input, "GLOBAL%PROJECT_NAME", c_val=project_name)
     224          54 :       CALL section_vals_val_get(input, "GLOBAL%OUTPUT_FILE_NAME", c_val=output_file_name)
     225             :       CALL rep_env_create(rep_env, para_env=para_env, input=input, &
     226          54 :                           input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=.TRUE.)
     227          54 :       CALL section_vals_val_set(input, "GLOBAL%PROJECT_NAME", c_val=TRIM(project_name))
     228          54 :       IF (LEN_TRIM(output_file_name) .GT. 0) THEN
     229           0 :          CALL section_vals_val_set(input, "GLOBAL%OUTPUT_FILE_NAME", c_val=TRIM(output_file_name))
     230             :       ELSE
     231          54 :          CALL section_vals_val_unset(input, "GLOBAL%OUTPUT_FILE_NAME")
     232             :       END IF
     233          54 :       IF (.NOT. ASSOCIATED(rep_env)) RETURN
     234             : 
     235          54 :       NULLIFY (pint_env%logger)
     236          54 :       pint_env%logger => cp_get_default_logger()
     237          54 :       CALL cp_add_iter_level(pint_env%logger%iter_info, "PINT")
     238             : 
     239          54 :       NULLIFY (pint_env%replicas, pint_env%input, pint_env%staging_env, &
     240          54 :                pint_env%normalmode_env, pint_env%propagator)
     241          54 :       pint_env%p = nrep
     242          54 :       pint_env%replicas => rep_env
     243          54 :       pint_env%ndim = rep_env%ndim
     244          54 :       pint_env%input => input
     245             : 
     246          54 :       CALL section_vals_retain(pint_env%input)
     247             : 
     248             :       ! get first step, last step, number of steps, etc
     249             :       CALL section_vals_val_get(input, "MOTION%PINT%ITERATION", &
     250          54 :                                 i_val=itmp)
     251          54 :       pint_env%first_step = itmp
     252             :       CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
     253          54 :                                 explicit=explicit)
     254          54 :       IF (explicit) THEN
     255             :          CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
     256           0 :                                    i_val=itmp)
     257           0 :          pint_env%last_step = itmp
     258           0 :          pint_env%num_steps = pint_env%last_step - pint_env%first_step
     259             :       ELSE
     260             :          CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", &
     261          54 :                                    i_val=itmp)
     262          54 :          pint_env%num_steps = itmp
     263          54 :          pint_env%last_step = pint_env%first_step + pint_env%num_steps
     264             :       END IF
     265             : 
     266             :       CALL section_vals_val_get(pint_section, "DT", &
     267          54 :                                 r_val=pint_env%dt)
     268          54 :       pint_env%t = pint_env%first_step*pint_env%dt
     269             : 
     270          54 :       CALL section_vals_val_get(pint_section, "nrespa", i_val=pint_env%nrespa)
     271          54 :       CALL section_vals_val_get(pint_section, "Temp", r_val=pint_env%kT)
     272             :       CALL section_vals_val_get(pint_section, "T_TOL", &
     273          54 :                                 r_val=pint_env%t_tol)
     274             : 
     275          54 :       CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator)
     276             : 
     277          54 :       ALLOCATE (pint_env%propagator)
     278             :       CALL section_vals_val_get(pint_section, "propagator", &
     279          54 :                                 i_val=pint_env%propagator%prop_kind)
     280             :       !Initialize simulation temperature depending on the propagator
     281             :       !As well as the scaling factor for the physical potential
     282          54 :       IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
     283          16 :          pint_env%propagator%temp_phys2sim = REAL(pint_env%p, dp)
     284          16 :          pint_env%propagator%physpotscale = 1.0_dp
     285             :       ELSE
     286          38 :          pint_env%propagator%temp_phys2sim = 1.0_dp
     287          38 :          pint_env%propagator%physpotscale = 1.0_dp/REAL(pint_env%p, dp)
     288             :       END IF
     289          54 :       pint_env%propagator%temp_sim2phys = 1.0_dp/pint_env%propagator%temp_phys2sim
     290          54 :       pint_env%kT = pint_env%kT*pint_env%propagator%temp_phys2sim
     291             : 
     292             :       CALL section_vals_val_get(pint_section, "transformation", &
     293          54 :                                 i_val=pint_env%transform)
     294             : 
     295          54 :       IF ((pint_env%propagator%prop_kind == propagator_cmd) .AND. &
     296             :           (pint_env%transform /= transformation_normal)) THEN
     297           0 :          CPABORT("CMD Propagator without normal modes not implemented!")
     298             :       END IF
     299             : 
     300          54 :       NULLIFY (pint_env%tx, pint_env%tv, pint_env%tv_t, pint_env%tv_old, pint_env%tv_new, pint_env%tf)
     301             : 
     302          54 :       pint_env%nnos = 0
     303          54 :       pint_env%pimd_thermostat = thermostat_none
     304          54 :       nose_section => section_vals_get_subs_vals(input, "MOTION%PINT%NOSE")
     305          54 :       CALL section_vals_get(nose_section, explicit=explicit)
     306          54 :       IF (explicit) THEN
     307          26 :          IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
     308           0 :             CPABORT("RPMD Propagator with Nose-thermostat not implemented!")
     309             :          END IF
     310          26 :          CALL section_vals_val_get(nose_section, "nnos", i_val=pint_env%nnos)
     311          26 :          IF (pint_env%nnos > 0) THEN
     312          26 :             pint_env%pimd_thermostat = thermostat_nose
     313             :             ALLOCATE ( &
     314             :                pint_env%tx(pint_env%nnos, pint_env%p, pint_env%ndim), &
     315             :                pint_env%tv(pint_env%nnos, pint_env%p, pint_env%ndim), &
     316             :                pint_env%tv_t(pint_env%nnos, pint_env%p, pint_env%ndim), &
     317             :                pint_env%tv_old(pint_env%nnos, pint_env%p, pint_env%ndim), &
     318             :                pint_env%tv_new(pint_env%nnos, pint_env%p, pint_env%ndim), &
     319         650 :                pint_env%tf(pint_env%nnos, pint_env%p, pint_env%ndim))
     320       88244 :             pint_env%tx = 0._dp
     321       88244 :             pint_env%tv = 0._dp
     322       88244 :             pint_env%tv_t = 0._dp
     323       88244 :             pint_env%tv_old = 0._dp
     324       88244 :             pint_env%tv_new = 0._dp
     325       88244 :             pint_env%tf = 0._dp
     326             :          END IF
     327             :       END IF
     328             : 
     329          54 :       pint_env%beta = 1._dp/(pint_env%kT*pint_env%propagator%temp_sim2phys)
     330             : !TODO
     331             : ! v_tol not in current input structure
     332             : ! should also probably be part of nose_section
     333             : !       CALL section_vals_val_get(transform_section,"v_tol_nose",r_val=pint_env%v_tol)
     334             : !MK ... but we have to initialise v_tol
     335          54 :       pint_env%v_tol = 0.0_dp ! to be fixed
     336             : 
     337             :       pint_env%randomG = rng_stream_type( &
     338             :                          name="pint_randomG", &
     339             :                          distribution_type=GAUSSIAN, &
     340          54 :                          extended_precision=.TRUE.)
     341             : 
     342         162 :       ALLOCATE (pint_env%e_pot_bead(pint_env%p))
     343         366 :       pint_env%e_pot_bead = 0._dp
     344          54 :       pint_env%e_pot_h = 0._dp
     345          54 :       pint_env%e_kin_beads = 0._dp
     346          54 :       pint_env%e_pot_t = 0._dp
     347          54 :       pint_env%e_gle = 0._dp
     348          54 :       pint_env%e_pile = 0._dp
     349          54 :       pint_env%e_piglet = 0._dp
     350          54 :       pint_env%e_qtb = 0._dp
     351          54 :       pint_env%e_kin_t = 0._dp
     352         270 :       pint_env%energy(:) = 0.0_dp
     353             : 
     354             : !TODO: rearrange to use standard nose hoover chain functions/data types
     355             : 
     356             :       ALLOCATE ( &
     357             :          pint_env%x(pint_env%p, pint_env%ndim), &
     358             :          pint_env%v(pint_env%p, pint_env%ndim), &
     359             :          pint_env%f(pint_env%p, pint_env%ndim), &
     360             :          pint_env%external_f(pint_env%p, pint_env%ndim), &
     361             :          pint_env%ux(pint_env%p, pint_env%ndim), &
     362             :          pint_env%ux_t(pint_env%p, pint_env%ndim), &
     363             :          pint_env%uv(pint_env%p, pint_env%ndim), &
     364             :          pint_env%uv_t(pint_env%p, pint_env%ndim), &
     365             :          pint_env%uv_new(pint_env%p, pint_env%ndim), &
     366             :          pint_env%uf(pint_env%p, pint_env%ndim), &
     367             :          pint_env%uf_h(pint_env%p, pint_env%ndim), &
     368             :          pint_env%centroid(pint_env%ndim), &
     369             :          pint_env%rtmp_ndim(pint_env%ndim), &
     370             :          pint_env%rtmp_natom(pint_env%ndim/3), &
     371        2160 :          STAT=stat)
     372          54 :       CPASSERT(stat == 0)
     373      362232 :       pint_env%x = 0._dp
     374      362232 :       pint_env%v = 0._dp
     375      362232 :       pint_env%f = 0._dp
     376      362232 :       pint_env%external_f = 0._dp
     377      362232 :       pint_env%ux = 0._dp
     378      362232 :       pint_env%ux_t = 0._dp
     379      362232 :       pint_env%uv = 0._dp
     380      362232 :       pint_env%uv_t = 0._dp
     381      362232 :       pint_env%uv_new = 0._dp
     382      362232 :       pint_env%uf = 0._dp
     383      362232 :       pint_env%uf_h = 0._dp
     384       64944 :       pint_env%centroid(:) = 0.0_dp
     385       64944 :       pint_env%rtmp_ndim = 0._dp
     386       21684 :       pint_env%rtmp_natom = 0._dp
     387          54 :       pint_env%time_per_step = 0.0_dp
     388             : 
     389          54 :       IF (pint_env%transform == transformation_stage) THEN
     390             :          transform_section => section_vals_get_subs_vals(input, &
     391           0 :                                                          "MOTION%PINT%STAGING")
     392           0 :          ALLOCATE (pint_env%staging_env)
     393             :          CALL staging_env_create(pint_env%staging_env, transform_section, &
     394           0 :                                  p=pint_env%p, kT=pint_env%kT)
     395             :       ELSE
     396             :          transform_section => section_vals_get_subs_vals(input, &
     397          54 :                                                          "MOTION%PINT%NORMALMODE")
     398          54 :          ALLOCATE (pint_env%normalmode_env)
     399             :          CALL normalmode_env_create(pint_env%normalmode_env, &
     400          54 :                                     transform_section, p=pint_env%p, kT=pint_env%kT, propagator=pint_env%propagator%prop_kind)
     401          54 :          IF (para_env%is_source()) THEN
     402          27 :             IF (pint_env%harm_integrator == integrate_numeric) THEN
     403         114 :                IF (10.0_dp*pint_env%dt/REAL(pint_env%nrespa, dp) > &
     404             :                    twopi/(pint_env%p*SQRT(MAXVAL(pint_env%normalmode_env%lambda))* &
     405             :                           pint_env%normalmode_env%modefactor)) THEN
     406             :                   msg = "PINT WARNING| Number of RESPA steps to small "// &
     407           0 :                         "to integrate the harmonic springs."
     408           0 :                   CPWARN(msg)
     409             :                END IF
     410             :             END IF
     411             :          END IF
     412             :       END IF
     413         162 :       ALLOCATE (pint_env%mass(pint_env%ndim))
     414             :       CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
     415          54 :                               f_env=f_env)
     416          54 :       CALL force_env_get(force_env=f_env%force_env, subsys=subsys)
     417          54 :       CALL cp_subsys_get(subsys, particles=particles, gci=gci)
     418             : 
     419             : !TODO length of pint_env%mass is redundant
     420          54 :       idim = 0
     421       21684 :       DO iat = 1, pint_env%ndim/3
     422       21630 :          CALL get_atomic_kind(particles%els(iat)%atomic_kind, mass=mass)
     423       86574 :          DO idir = 1, 3
     424       64890 :             idim = idim + 1
     425       86520 :             pint_env%mass(idim) = mass
     426             :          END DO
     427             :       END DO
     428          54 :       CALL f_env_rm_defaults(f_env, ierr)
     429          54 :       CPASSERT(ierr == 0)
     430             : 
     431             :       ALLOCATE (pint_env%Q(pint_env%p), &
     432             :                 pint_env%mass_beads(pint_env%p, pint_env%ndim), &
     433         486 :                 pint_env%mass_fict(pint_env%p, pint_env%ndim))
     434          54 :       IF (pint_env%transform == transformation_stage) THEN
     435             :          CALL staging_init_masses(pint_env%staging_env, mass=pint_env%mass, &
     436             :                                   mass_beads=pint_env%mass_beads, mass_fict=pint_env%mass_fict, &
     437           0 :                                   Q=pint_env%Q)
     438             :       ELSE
     439             :          CALL normalmode_init_masses(pint_env%normalmode_env, &
     440             :                                      mass=pint_env%mass, mass_beads=pint_env%mass_beads, &
     441          54 :                                      mass_fict=pint_env%mass_fict, Q=pint_env%Q)
     442             :       END IF
     443             : 
     444          54 :       NULLIFY (pint_env%gle)
     445          54 :       gle_section => section_vals_get_subs_vals(input, "MOTION%PINT%GLE")
     446          54 :       CALL section_vals_get(gle_section, explicit=explicit)
     447          54 :       IF (explicit) THEN
     448           2 :          ALLOCATE (pint_env%gle)
     449             :          CALL gle_init(pint_env%gle, dt=pint_env%dt/pint_env%nrespa, temp=pint_env%kT, &
     450           2 :                        section=gle_section)
     451           2 :          IF (pint_env%pimd_thermostat == thermostat_none .AND. pint_env%gle%ndim .GT. 0) THEN
     452           2 :             pint_env%pimd_thermostat = thermostat_gle
     453             : 
     454             :             ! initialize a GLE with ALL degrees of freedom on node 0,
     455             :             ! as it seems to me that here everything but force eval is replicated
     456           2 :             pint_env%gle%loc_num_gle = pint_env%p*pint_env%ndim
     457           2 :             pint_env%gle%glob_num_gle = pint_env%gle%loc_num_gle
     458           6 :             ALLOCATE (pint_env%gle%map_info%index(pint_env%gle%loc_num_gle))
     459           2 :             CPASSERT(stat == 0)
     460       18434 :             DO itmp = 1, pint_env%gle%loc_num_gle
     461       18434 :                pint_env%gle%map_info%index(itmp) = itmp
     462             :             END DO
     463           2 :             CALL gle_thermo_create(pint_env%gle, pint_env%gle%loc_num_gle)
     464             : 
     465             :             ! here we should have read a_mat and c_mat;
     466             :             !we can therefore compute the matrices needed for the propagator
     467             :             ! deterministic part of the propagator
     468             :             CALL gle_matrix_exp((-pint_env%dt/pint_env%nrespa*0.5_dp)*pint_env%gle%a_mat, &
     469          62 :                                 pint_env%gle%ndim, 15, 15, pint_env%gle%gle_t)
     470             :             ! stochastic part
     471           8 :             CALL gle_cholesky_stab(pint_env%gle%c_mat - MATMUL(pint_env%gle%gle_t, &
     472           8 :                                                                MATMUL(pint_env%gle%c_mat, TRANSPOSE(pint_env%gle%gle_t))), &
     473        2304 :                                    pint_env%gle%gle_s, pint_env%gle%ndim)
     474             :             ! and initialize the additional momenta
     475           2 :             CALL pint_gle_init(pint_env)
     476             :          END IF
     477             :       END IF
     478             : 
     479             :       !Setup pile thermostat
     480          54 :       NULLIFY (pint_env%pile_therm)
     481          54 :       pile_section => section_vals_get_subs_vals(input, "MOTION%PINT%PILE")
     482          54 :       CALL section_vals_get(pile_section, explicit=explicit)
     483          54 :       IF (explicit) THEN
     484           8 :          CALL cite_reference(Ceriotti2010)
     485             :          CALL section_vals_val_get(pint_env%input, &
     486             :                                    "MOTION%PINT%INIT%THERMOSTAT_SEED", &
     487           8 :                                    i_val=pint_env%thermostat_rng_seed)
     488           8 :          IF (pint_env%pimd_thermostat == thermostat_none) THEN
     489           8 :             pint_env%pimd_thermostat = thermostat_pile
     490         200 :             ALLOCATE (pint_env%pile_therm)
     491             :             CALL pint_pile_init(pile_therm=pint_env%pile_therm, &
     492             :                                 pint_env=pint_env, &
     493             :                                 normalmode_env=pint_env%normalmode_env, &
     494           8 :                                 section=pile_section)
     495             :          ELSE
     496           0 :             CPABORT("PILE thermostat can't be used with another thermostat.")
     497             :          END IF
     498             :       END IF
     499             : 
     500             :       !Setup PIGLET thermostat
     501          54 :       NULLIFY (pint_env%piglet_therm)
     502          54 :       piglet_section => section_vals_get_subs_vals(input, "MOTION%PINT%PIGLET")
     503          54 :       CALL section_vals_get(piglet_section, explicit=explicit)
     504          54 :       IF (explicit) THEN
     505           2 :          CALL cite_reference(Ceriotti2012)
     506             :          CALL section_vals_val_get(pint_env%input, &
     507             :                                    "MOTION%PINT%INIT%THERMOSTAT_SEED", &
     508           2 :                                    i_val=pint_env%thermostat_rng_seed)
     509           2 :          IF (pint_env%pimd_thermostat == thermostat_none) THEN
     510           2 :             pint_env%pimd_thermostat = thermostat_piglet
     511          50 :             ALLOCATE (pint_env%piglet_therm)
     512             :             CALL pint_piglet_create(pint_env%piglet_therm, &
     513             :                                     pint_env, &
     514           2 :                                     piglet_section)
     515             :             CALL pint_piglet_init(pint_env%piglet_therm, &
     516             :                                   pint_env, &
     517             :                                   piglet_section, &
     518           2 :                                   dt=pint_env%dt, para_env=para_env)
     519             :          ELSE
     520           0 :             CPABORT("PILE thermostat can't be used with another thermostat.")
     521             :          END IF
     522             :       END IF
     523             : 
     524             :       !Setup qtb thermostat
     525          54 :       NULLIFY (pint_env%qtb_therm)
     526          54 :       qtb_section => section_vals_get_subs_vals(input, "MOTION%PINT%QTB")
     527          54 :       CALL section_vals_get(qtb_section, explicit=explicit)
     528          54 :       IF (explicit) THEN
     529           6 :          CALL cite_reference(Brieuc2016)
     530             :          CALL section_vals_val_get(pint_env%input, &
     531             :                                    "MOTION%PINT%INIT%THERMOSTAT_SEED", &
     532           6 :                                    i_val=pint_env%thermostat_rng_seed)
     533           6 :          IF (pint_env%pimd_thermostat == thermostat_none) THEN
     534           6 :             pint_env%pimd_thermostat = thermostat_qtb
     535             :             CALL pint_qtb_init(qtb_therm=pint_env%qtb_therm, &
     536             :                                pint_env=pint_env, &
     537             :                                normalmode_env=pint_env%normalmode_env, &
     538           6 :                                section=qtb_section)
     539             :          ELSE
     540           0 :             CPABORT("QTB thermostat can't be used with another thermostat.")
     541             :          END IF
     542             :       END IF
     543             : 
     544             :       !Initialize integrator scheme
     545          54 :       CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator)
     546          54 :       IF (pint_env%harm_integrator == integrate_exact) THEN
     547          20 :          IF (pint_env%pimd_thermostat == thermostat_nose) THEN
     548             :             WRITE (UNIT=msg, FMT=*) "PINT WARNING| Nose Thermostat only available in "// &
     549           0 :                "the numeric harmonic integrator. Switching to numeric harmonic integrator."
     550           0 :             CPWARN(msg)
     551           0 :             pint_env%harm_integrator = integrate_numeric
     552             :          END IF
     553          20 :          IF (pint_env%pimd_thermostat == thermostat_gle) THEN
     554             :             WRITE (UNIT=msg, FMT=*) "PINT WARNING| GLE Thermostat only available in "// &
     555           0 :                "the numeric harmonic integrator. Switching to numeric harmonic integrator."
     556           0 :             CPWARN(msg)
     557           0 :             pint_env%harm_integrator = integrate_numeric
     558             :          END IF
     559          34 :       ELSE IF (pint_env%harm_integrator == integrate_numeric) THEN
     560          34 :          IF (pint_env%pimd_thermostat == thermostat_pile) THEN
     561             :             WRITE (UNIT=msg, FMT=*) "PINT WARNING| PILE Thermostat only available in "// &
     562           2 :                "the exact harmonic integrator. Switching to exact harmonic integrator."
     563           2 :             CPWARN(msg)
     564           2 :             pint_env%harm_integrator = integrate_exact
     565             :          END IF
     566          34 :          IF (pint_env%pimd_thermostat == thermostat_piglet) THEN
     567             :             WRITE (UNIT=msg, FMT=*) "PINT WARNING| PIGLET Thermostat only available in "// &
     568           0 :                "the exact harmonic integrator. Switching to exact harmonic integrator."
     569           0 :             CPWARN(msg)
     570           0 :             pint_env%harm_integrator = integrate_exact
     571             :          END IF
     572          34 :          IF (pint_env%pimd_thermostat == thermostat_qtb) THEN
     573             :             WRITE (UNIT=msg, FMT=*) "PINT WARNING| QTB Thermostat only available in "// &
     574           0 :                "the exact harmonic integrator. Switching to exact harmonic integrator."
     575           0 :             CPWARN(msg)
     576           0 :             pint_env%harm_integrator = integrate_exact
     577             :          END IF
     578             :       END IF
     579             : 
     580          54 :       IF (pint_env%harm_integrator == integrate_exact) THEN
     581          22 :          IF (pint_env%nrespa /= 1) THEN
     582          18 :             pint_env%nrespa = 1
     583          18 :             WRITE (UNIT=msg, FMT=*) "PINT WARNING| Adjusting NRESPA to 1 for exact harmonic integration."
     584          18 :             CPWARN(msg)
     585             :          END IF
     586          22 :          NULLIFY (pint_env%wsinex)
     587          66 :          ALLOCATE (pint_env%wsinex(pint_env%p))
     588          22 :          NULLIFY (pint_env%iwsinex)
     589          66 :          ALLOCATE (pint_env%iwsinex(pint_env%p))
     590          22 :          NULLIFY (pint_env%cosex)
     591          66 :          ALLOCATE (pint_env%cosex(pint_env%p))
     592          22 :          dt = pint_env%dt/REAL(pint_env%nrespa, KIND=dp)
     593             :          !Centroid mode shoud not be propagated
     594          22 :          pint_env%wsinex(1) = 0.0_dp
     595          22 :          pint_env%iwsinex(1) = dt
     596          22 :          pint_env%cosex(1) = 1.0_dp
     597         172 :          DO ibead = 2, pint_env%p
     598         150 :             omega = SQRT(pint_env%normalmode_env%lambda(ibead))
     599         150 :             pint_env%wsinex(ibead) = SIN(omega*dt)*omega
     600         150 :             pint_env%iwsinex(ibead) = SIN(omega*dt)/omega
     601         172 :             pint_env%cosex(ibead) = COS(omega*dt)
     602             :          END DO
     603             :       END IF
     604             : 
     605             :       CALL section_vals_val_get(pint_section, "FIX_CENTROID_POS", &
     606          54 :                                 l_val=ltmp)
     607          54 :       IF (ltmp .AND. (pint_env%transform .EQ. transformation_normal)) THEN
     608           0 :          pint_env%first_propagated_mode = 2
     609             :       ELSE
     610          54 :          pint_env%first_propagated_mode = 1
     611             :       END IF
     612             : 
     613             :       ! Constraint information:
     614          54 :       NULLIFY (pint_env%simpar)
     615             :       constraint_section => section_vals_get_subs_vals(pint_env%input, &
     616          54 :                                                        "MOTION%CONSTRAINT")
     617          54 :       CALL section_vals_get(constraint_section, explicit=explicit)
     618          54 :       CALL create_simpar_type(pint_env%simpar)
     619          54 :       pint_env%simpar%constraint = explicit
     620          54 :       pint_env%kTcorr = 1.0_dp
     621             : 
     622             :       ! Determine if beadwise constraints are activated
     623          54 :       pint_env%beadwise_constraints = .FALSE.
     624             :       CALL section_vals_val_get(constraint_section, "PIMD_BEADWISE_CONSTRAINT", &
     625          54 :                                 l_val=pint_env%beadwise_constraints)
     626          54 :       IF (pint_env%beadwise_constraints) THEN
     627           2 :          CALL pint_write_line("Using beadwise constraints")
     628             :       ELSE
     629          52 :          CALL pint_write_line("Using centroid constraints")
     630             :       END IF
     631             : 
     632          54 :       IF (explicit) THEN
     633             :          ! Staging not supported yet, since lowest mode is assumed
     634             :          ! to be related to centroid
     635           6 :          IF (pint_env%transform == transformation_stage) THEN
     636           0 :             CPABORT("Constraints are not supported for staging transformation")
     637             :          END IF
     638             : 
     639             :          ! Check thermostats that are not supported:
     640           6 :          IF (pint_env%pimd_thermostat == thermostat_gle) THEN
     641             :             WRITE (UNIT=msg, FMT=*) "GLE Thermostat not supported for "// &
     642           0 :                "constraints. Switch to NOSE for numeric integration."
     643           0 :             CPABORT(msg)
     644             :          END IF
     645             :          ! Warn for NOSE
     646           6 :          IF (pint_env%pimd_thermostat == thermostat_nose) THEN
     647             :             !Beadwise constraints not supported
     648           2 :             IF (pint_env%beadwise_constraints) THEN
     649           0 :                CPABORT("Beadwise constraints are not supported for NOSE Thermostat.")
     650             :                !Centroid constraints supported
     651             :             ELSE
     652             :                WRITE (UNIT=msg, FMT=*) "PINT WARNING| Nose Thermostat set to "// &
     653           2 :                   "zero for constrained atoms. Careful interpretation of temperature."
     654           2 :                CPWARN(msg)
     655             :                WRITE (UNIT=msg, FMT=*) "PINT WARNING| Lagrange multipliers are "// &
     656           2 :                   "are printed every RESPA step and need to be treated carefully."
     657           2 :                CPWARN(msg)
     658             :             END IF
     659             :          END IF
     660             : 
     661             :          CALL section_vals_val_get(constraint_section, "SHAKE_TOLERANCE", &
     662           6 :                                    r_val=pint_env%simpar%shake_tol)
     663             :          pint_env%simpar%info_constraint = cp_print_key_unit_nr(pint_env%logger, &
     664             :                                                                 constraint_section, &
     665             :                                                                 "CONSTRAINT_INFO", &
     666             :                                                                 extension=".shakeLog", &
     667           6 :                                                                 log_filename=.FALSE.)
     668             :          pint_env%simpar%lagrange_multipliers = cp_print_key_unit_nr(pint_env%logger, &
     669             :                                                                      constraint_section, &
     670             :                                                                      "LAGRANGE_MULTIPLIERS", &
     671             :                                                                      extension=".LagrangeMultLog", &
     672           6 :                                                                      log_filename=.FALSE.)
     673             :          pint_env%simpar%dump_lm = &
     674             :             BTEST(cp_print_key_should_output(pint_env%logger%iter_info, &
     675             :                                              constraint_section, &
     676           6 :                                              "LAGRANGE_MULTIPLIERS"), cp_p_file)
     677             : 
     678             :          ! Determine constrained atoms:
     679           6 :          pint_env%n_atoms_constraints = 0
     680          12 :          DO ig = 1, gci%ncolv%ntot
     681             :             ! Double counts, if the same atom is involved in different collective variables
     682          12 :             pint_env%n_atoms_constraints = pint_env%n_atoms_constraints + SIZE(gci%colv_list(ig)%i_atoms)
     683             :          END DO
     684             : 
     685          18 :          ALLOCATE (pint_env%atoms_constraints(pint_env%n_atoms_constraints))
     686           6 :          icont = 0
     687          12 :          DO ig = 1, gci%ncolv%ntot
     688          24 :             DO iat = 1, SIZE(gci%colv_list(ig)%i_atoms)
     689          12 :                icont = icont + 1
     690          18 :                pint_env%atoms_constraints(icont) = gci%colv_list(ig)%i_atoms(iat)
     691             :             END DO
     692             :          END DO
     693             : 
     694             :          ! Set the correction to the temperature due to the frozen degrees of freedom in NOSE:
     695             :          CALL section_vals_val_get(pint_section, "kT_CORRECTION", &
     696           6 :                                    l_val=ltmp)
     697           6 :          IF (ltmp) THEN
     698           0 :             pint_env%kTcorr = 1.0_dp + REAL(3*pint_env%n_atoms_constraints, dp)/(REAL(pint_env%ndim, dp)*REAL(pint_env%p, dp))
     699             :          END IF
     700             :       END IF
     701             : 
     702          54 :       CALL timestop(handle)
     703             : 
     704         594 :    END SUBROUTINE pint_create
     705             : 
     706             : ! ***************************************************************************
     707             : !> \brief Release a path integral environment
     708             : !> \param pint_env the pint_env to release
     709             : !> \par History
     710             : !>      Added normal mode transformation [hforbert]
     711             : !> \author Fawzi Mohamed
     712             : ! **************************************************************************************************
     713          54 :    SUBROUTINE pint_release(pint_env)
     714             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     715             : 
     716          54 :       CALL rep_env_release(pint_env%replicas)
     717          54 :       CALL section_vals_release(pint_env%input)
     718          54 :       IF (ASSOCIATED(pint_env%staging_env)) THEN
     719           0 :          CALL staging_release(pint_env%staging_env)
     720           0 :          DEALLOCATE (pint_env%staging_env)
     721             :       END IF
     722          54 :       IF (ASSOCIATED(pint_env%normalmode_env)) THEN
     723          54 :          CALL normalmode_release(pint_env%normalmode_env)
     724          54 :          DEALLOCATE (pint_env%normalmode_env)
     725             :       END IF
     726             : 
     727          54 :       DEALLOCATE (pint_env%mass)
     728          54 :       DEALLOCATE (pint_env%e_pot_bead)
     729             : 
     730          54 :       DEALLOCATE (pint_env%x)
     731          54 :       DEALLOCATE (pint_env%v)
     732          54 :       DEALLOCATE (pint_env%f)
     733          54 :       DEALLOCATE (pint_env%external_f)
     734          54 :       DEALLOCATE (pint_env%mass_beads)
     735          54 :       DEALLOCATE (pint_env%mass_fict)
     736          54 :       DEALLOCATE (pint_env%ux)
     737          54 :       DEALLOCATE (pint_env%ux_t)
     738          54 :       DEALLOCATE (pint_env%uv)
     739          54 :       DEALLOCATE (pint_env%uv_t)
     740          54 :       DEALLOCATE (pint_env%uv_new)
     741          54 :       DEALLOCATE (pint_env%uf)
     742          54 :       DEALLOCATE (pint_env%uf_h)
     743          54 :       DEALLOCATE (pint_env%centroid)
     744          54 :       DEALLOCATE (pint_env%rtmp_ndim)
     745          54 :       DEALLOCATE (pint_env%rtmp_natom)
     746          54 :       DEALLOCATE (pint_env%propagator)
     747             : 
     748          54 :       IF (pint_env%simpar%constraint) THEN
     749           6 :          DEALLOCATE (pint_env%atoms_constraints)
     750             :       END IF
     751          54 :       CALL release_simpar_type(pint_env%simpar)
     752             : 
     753          54 :       IF (pint_env%harm_integrator == integrate_exact) THEN
     754          22 :          DEALLOCATE (pint_env%wsinex)
     755          22 :          DEALLOCATE (pint_env%iwsinex)
     756          22 :          DEALLOCATE (pint_env%cosex)
     757             :       END IF
     758             : 
     759          80 :       SELECT CASE (pint_env%pimd_thermostat)
     760             :       CASE (thermostat_nose)
     761          26 :          DEALLOCATE (pint_env%tx)
     762          26 :          DEALLOCATE (pint_env%tv)
     763          26 :          DEALLOCATE (pint_env%tv_t)
     764          26 :          DEALLOCATE (pint_env%tv_old)
     765          26 :          DEALLOCATE (pint_env%tv_new)
     766          26 :          DEALLOCATE (pint_env%tf)
     767             :       CASE (thermostat_gle)
     768           2 :          CALL gle_dealloc(pint_env%gle)
     769             :       CASE (thermostat_pile)
     770           8 :          CALL pint_pile_release(pint_env%pile_therm)
     771           8 :          DEALLOCATE (pint_env%pile_therm)
     772             :       CASE (thermostat_piglet)
     773           2 :          CALL pint_piglet_release(pint_env%piglet_therm)
     774           2 :          DEALLOCATE (pint_env%piglet_therm)
     775             :       CASE (thermostat_qtb)
     776           6 :          CALL pint_qtb_release(pint_env%qtb_therm)
     777          60 :          DEALLOCATE (pint_env%qtb_therm)
     778             :       END SELECT
     779             : 
     780          54 :       DEALLOCATE (pint_env%Q)
     781             : 
     782          54 :    END SUBROUTINE pint_release
     783             : 
     784             : ! ***************************************************************************
     785             : !> \brief Tests the path integral methods
     786             : !> \param para_env parallel environment
     787             : !> \param input the input to test
     788             : !> \param input_declaration ...
     789             : !> \author fawzi
     790             : ! **************************************************************************************************
     791           0 :    SUBROUTINE pint_test(para_env, input, input_declaration)
     792             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     793             :       TYPE(section_vals_type), POINTER                   :: input
     794             :       TYPE(section_type), POINTER                        :: input_declaration
     795             : 
     796             :       INTEGER                                            :: i, ib, idim, unit_nr
     797             :       REAL(kind=dp)                                      :: c, e_h, err
     798           0 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: x1
     799             :       TYPE(pint_env_type)                                :: pint_env
     800             : 
     801           0 :       CPASSERT(ASSOCIATED(para_env))
     802           0 :       CPASSERT(ASSOCIATED(input))
     803           0 :       CPASSERT(para_env%is_valid())
     804           0 :       CPASSERT(input%ref_count > 0)
     805           0 :       unit_nr = cp_logger_get_default_io_unit()
     806           0 :       CALL pint_create(pint_env, input, input_declaration, para_env)
     807           0 :       ALLOCATE (x1(pint_env%ndim, pint_env%p))
     808           0 :       x1(:, :) = pint_env%x
     809           0 :       CALL pint_x2u(pint_env)
     810           0 :       pint_env%x = 0._dp
     811           0 :       CALL pint_u2x(pint_env)
     812           0 :       err = 0._dp
     813           0 :       DO i = 1, pint_env%ndim
     814           0 :          err = MAX(err, ABS(x1(1, i) - pint_env%x(1, i)))
     815             :       END DO
     816           0 :       IF (unit_nr > 0) WRITE (unit_nr, *) "diff_r1="//cp_to_string(err)
     817             : 
     818           0 :       CALL pint_calc_uf_h(pint_env, e_h=e_h)
     819           0 :       c = -pint_env%staging_env%w_p**2
     820           0 :       pint_env%f = 0._dp
     821           0 :       DO idim = 1, pint_env%ndim
     822           0 :          DO ib = 1, pint_env%p
     823             :             pint_env%f(ib, idim) = pint_env%f(ib, idim) + &
     824             :                                    c*(2._dp*pint_env%x(ib, idim) &
     825             :                                       - pint_env%x(MODULO(ib - 2, pint_env%p) + 1, idim) &
     826           0 :                                       - pint_env%x(MODULO(ib, pint_env%p) + 1, idim))
     827             :          END DO
     828             :       END DO
     829           0 :       CALL pint_f2uf(pint_env)
     830           0 :       err = 0._dp
     831           0 :       DO idim = 1, pint_env%ndim
     832           0 :          DO ib = 1, pint_env%p
     833           0 :             err = MAX(err, ABS(pint_env%uf(ib, idim) - pint_env%uf_h(ib, idim)))
     834             :          END DO
     835             :       END DO
     836           0 :       IF (unit_nr > 0) WRITE (unit_nr, *) "diff_f_h="//cp_to_string(err)
     837             : 
     838           0 :    END SUBROUTINE pint_test
     839             : 
     840             : ! ***************************************************************************
     841             : !> \brief  Perform a path integral simulation
     842             : !> \param  para_env parallel environment
     843             : !> \param  input the input to test
     844             : !> \param input_declaration ...
     845             : !> \param globenv ...
     846             : !> \par    History
     847             : !>         2003-11 created [fawzi]
     848             : !>         2009-12-14 globenv parameter added to handle soft exit
     849             : !>           requests [lwalewski]
     850             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     851             : !> \author Fawzi Mohamed
     852             : ! **************************************************************************************************
     853         192 :    SUBROUTINE do_pint_run(para_env, input, input_declaration, globenv)
     854             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     855             :       TYPE(section_vals_type), POINTER                   :: input
     856             :       TYPE(section_type), POINTER                        :: input_declaration
     857             :       TYPE(global_environment_type), POINTER             :: globenv
     858             : 
     859             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'do_pint_run'
     860             :       INTEGER, PARAMETER                                 :: helium_only_mid = 1, &
     861             :                                                             int_pot_scan_mid = 4, &
     862             :                                                             solute_only_mid = 2, &
     863             :                                                             solute_with_helium_mid = 3
     864             : 
     865             :       CHARACTER(len=default_string_length)               :: stmp
     866             :       INTEGER                                            :: handle, mode
     867             :       LOGICAL                                            :: explicit, helium_only, int_pot_scan, &
     868             :                                                             solvent_present
     869          64 :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     870             :       TYPE(pint_env_type)                                :: pint_env
     871             :       TYPE(section_vals_type), POINTER                   :: helium_section
     872             : 
     873          64 :       CALL timeset(routineN, handle)
     874             : 
     875          64 :       CPASSERT(ASSOCIATED(para_env))
     876          64 :       CPASSERT(ASSOCIATED(input))
     877          64 :       CPASSERT(para_env%is_valid())
     878          64 :       CPASSERT(input%ref_count > 0)
     879             : 
     880             :       ! check if helium solvent is present
     881          64 :       NULLIFY (helium_section)
     882             :       helium_section => section_vals_get_subs_vals(input, &
     883          64 :                                                    "MOTION%PINT%HELIUM")
     884          64 :       CALL section_vals_get(helium_section, explicit=explicit)
     885          64 :       IF (explicit) THEN
     886             :          CALL section_vals_val_get(helium_section, "_SECTION_PARAMETERS_", &
     887          24 :                                    l_val=solvent_present)
     888             :       ELSE
     889          40 :          solvent_present = .FALSE.
     890             :       END IF
     891             : 
     892             :       ! check if there is anything but helium
     893          64 :       IF (solvent_present) THEN
     894             :          CALL section_vals_val_get(helium_section, "HELIUM_ONLY", &
     895          24 :                                    l_val=helium_only)
     896             :       ELSE
     897          40 :          helium_only = .FALSE.
     898             :       END IF
     899             : 
     900             :       ! check wheather to perform solute-helium interaction pot scan
     901          64 :       IF (solvent_present) THEN
     902             :          CALL section_vals_val_get(helium_section, "INTERACTION_POT_SCAN", &
     903          24 :                                    l_val=int_pot_scan)
     904             :       ELSE
     905          40 :          int_pot_scan = .FALSE.
     906             :       END IF
     907             : 
     908             :       ! input consistency check
     909          64 :       IF (helium_only .AND. int_pot_scan) THEN
     910           0 :          stmp = "Options HELIUM_ONLY and INTERACTION_POT_SCAN are exclusive"
     911           0 :          CPABORT(stmp)
     912             :       END IF
     913             : 
     914             :       ! select mode of operation
     915             :       mode = 0
     916          64 :       IF (solvent_present) THEN
     917          24 :          IF (helium_only) THEN
     918          10 :             mode = helium_only_mid
     919             :          ELSE
     920          14 :             IF (int_pot_scan) THEN
     921           0 :                mode = int_pot_scan_mid
     922             :             ELSE
     923          14 :                mode = solute_with_helium_mid
     924             :             END IF
     925             :          END IF
     926             :       ELSE
     927          40 :          mode = solute_only_mid
     928             :       END IF
     929             : 
     930             :       ! perform the simulation according to the chosen mode
     931          10 :       SELECT CASE (mode)
     932             : 
     933             :       CASE (helium_only_mid)
     934          10 :          CALL helium_create(helium_env, input)
     935          10 :          CALL helium_init(helium_env, pint_env)
     936          10 :          CALL helium_do_run(helium_env, globenv)
     937          10 :          CALL helium_release(helium_env)
     938             : 
     939             :       CASE (solute_only_mid)
     940          40 :          CALL pint_create(pint_env, input, input_declaration, para_env)
     941          40 :          CALL pint_init(pint_env)
     942          40 :          CALL pint_do_run(pint_env, globenv)
     943          40 :          CALL pint_release(pint_env)
     944             : 
     945             :       CASE (int_pot_scan_mid)
     946           0 :          CALL pint_create(pint_env, input, input_declaration, para_env)
     947             : ! TODO only initialization of positions is necessary, but rep_env_calc_e_f called
     948             : ! from within pint_init_f does something to the replica environments which can not be
     949             : ! avoided (has something to do with f_env_add_defaults) so leaving for now..
     950           0 :          CALL pint_init(pint_env)
     951           0 :          CALL helium_create(helium_env, input, solute=pint_env)
     952           0 :          CALL pint_run_scan(pint_env, helium_env)
     953           0 :          CALL helium_release(helium_env)
     954           0 :          CALL pint_release(pint_env)
     955             : 
     956             :       CASE (solute_with_helium_mid)
     957          14 :          CALL pint_create(pint_env, input, input_declaration, para_env)
     958             :          ! init pint without helium forces (they are not yet initialized)
     959          14 :          CALL pint_init(pint_env)
     960             :          ! init helium with solute's positions (they are already initialized)
     961          14 :          CALL helium_create(helium_env, input, solute=pint_env)
     962          14 :          CALL helium_init(helium_env, pint_env)
     963             :          ! reinit pint forces with helium forces (they are now initialized)
     964          14 :          CALL pint_init_f(pint_env, helium_env=helium_env)
     965             : 
     966          14 :          CALL pint_do_run(pint_env, globenv, helium_env=helium_env)
     967          14 :          CALL helium_release(helium_env)
     968          14 :          CALL pint_release(pint_env)
     969             : 
     970             :       CASE DEFAULT
     971          64 :          CPABORT("Unknown mode ("//TRIM(ADJUSTL(cp_to_string(mode)))//")")
     972             :       END SELECT
     973             : 
     974          64 :       CALL timestop(handle)
     975             : 
     976        1600 :    END SUBROUTINE do_pint_run
     977             : 
     978             : ! ***************************************************************************
     979             : !> \brief  Reads the restart, initializes the beads, etc.
     980             : !> \param pint_env ...
     981             : !> \par    History
     982             : !>           11.2003 created [fawzi]
     983             : !>           actually ASSIGN input pointer [hforbert]
     984             : !>           2010-12-16 turned into a wrapper routine [lwalewski]
     985             : !> \author Fawzi Mohamed
     986             : ! **************************************************************************************************
     987          54 :    SUBROUTINE pint_init(pint_env)
     988             : 
     989             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     990             : 
     991          54 :       CALL pint_init_x(pint_env)
     992          54 :       CALL pint_init_v(pint_env)
     993          54 :       CALL pint_init_t(pint_env)
     994          54 :       CALL pint_init_f(pint_env)
     995             : 
     996          54 :    END SUBROUTINE pint_init
     997             : 
     998             : ! ***************************************************************************
     999             : !> \brief  Assign initial postions to the beads.
    1000             : !> \param pint_env ...
    1001             : !> \date   2010-12-15
    1002             : !> \author Lukasz Walewski
    1003             : !> \note  Initialization is done in the following way:
    1004             : !>           1. assign all beads with the same classical positions from
    1005             : !>              FORCE_EVAL (hot start)
    1006             : !>           2. spread the beads around classical positions as if they were
    1007             : !>              free particles (if requested)
    1008             : !>           3. replace positions generated in steps 1-2 with the explicit
    1009             : !>              ones if they are explicitly given in the input structure
    1010             : !>           4. apply Gaussian noise to the positions generated so far (if
    1011             : !>              requested)
    1012             : ! **************************************************************************************************
    1013          54 :    SUBROUTINE pint_init_x(pint_env)
    1014             : 
    1015             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1016             : 
    1017             :       CHARACTER(len=5*default_string_length)             :: msg, tmp
    1018             :       INTEGER                                            :: ia, ib, ic, idim, input_seed, n_rep_val
    1019             :       LOGICAL                                            :: done_init, done_levy, done_rand, &
    1020             :                                                             explicit, levycorr, ltmp
    1021             :       REAL(kind=dp)                                      :: tcorr, var
    1022             :       REAL(kind=dp), DIMENSION(3)                        :: x0
    1023             :       REAL(kind=dp), DIMENSION(3, 2)                     :: seed
    1024          54 :       REAL(kind=dp), DIMENSION(:), POINTER               :: bx, r_vals
    1025          54 :       TYPE(rng_stream_type), ALLOCATABLE                 :: rng_gaussian
    1026             :       TYPE(section_vals_type), POINTER                   :: input_section
    1027             : 
    1028       64944 :       DO idim = 1, pint_env%ndim
    1029      362232 :          DO ib = 1, pint_env%p
    1030      362178 :             pint_env%x(ib, idim) = pint_env%replicas%r(idim, ib)
    1031             :          END DO
    1032             :       END DO
    1033             : 
    1034          54 :       done_levy = .FALSE.
    1035             :       CALL section_vals_val_get(pint_env%input, &
    1036             :                                 "MOTION%PINT%INIT%LEVY_POS_SAMPLE", &
    1037          54 :                                 l_val=ltmp)
    1038             :       CALL section_vals_val_get(pint_env%input, &
    1039             :                                 "MOTION%PINT%INIT%LEVY_TEMP_FACTOR", &
    1040          54 :                                 r_val=tcorr)
    1041          54 :       IF (ltmp) THEN
    1042             : 
    1043           0 :          IF (pint_env%beadwise_constraints) THEN
    1044             :             WRITE (UNIT=msg, FMT=*) "Beadwise constraints are not supported for "// &
    1045             :                "the initialization of the beads as free particles. "// &
    1046           0 :                "Please use hot start (default)."
    1047           0 :             CPABORT(msg)
    1048             :          END IF
    1049             : 
    1050           0 :          NULLIFY (bx)
    1051           0 :          ALLOCATE (bx(3*pint_env%p))
    1052             :          CALL section_vals_val_get(pint_env%input, &
    1053           0 :                                    "MOTION%PINT%INIT%LEVY_SEED", i_val=input_seed)
    1054           0 :          seed(:, :) = REAL(input_seed, KIND=dp)
    1055             : !      seed(:,:) = next_rng_seed()
    1056             :          rng_gaussian = rng_stream_type( &
    1057             :                         name="tmp_rng_gaussian", &
    1058             :                         distribution_type=GAUSSIAN, &
    1059             :                         extended_precision=.TRUE., &
    1060           0 :                         seed=seed)
    1061             : 
    1062             :          CALL section_vals_val_get(pint_env%input, &
    1063             :                                    "MOTION%PINT%INIT%LEVY_CORRELATED", &
    1064           0 :                                    l_val=levycorr)
    1065             : 
    1066           0 :          IF (levycorr) THEN
    1067             : 
    1068             :             ! correlated Levy walk - the same path for all atoms
    1069           0 :             x0 = (/0.0_dp, 0.0_dp, 0.0_dp/)
    1070           0 :             CALL pint_levy_walk(x0, pint_env%p, 1.0_dp, bx, rng_gaussian)
    1071           0 :             idim = 0
    1072           0 :             DO ia = 1, pint_env%ndim/3
    1073           0 :                var = SQRT(1.0_dp/(pint_env%kT*tcorr*pint_env%mass(3*ia)))
    1074           0 :                DO ic = 1, 3
    1075           0 :                   idim = idim + 1
    1076           0 :                   DO ib = 1, pint_env%p
    1077           0 :                      pint_env%x(ib, idim) = pint_env%x(ib, idim) + bx(3*(ib - 1) + ic)*var
    1078             :                   END DO
    1079             :                END DO
    1080             :             END DO
    1081             : 
    1082             :          ELSE
    1083             : 
    1084             :             ! uncorrelated bead initialization - distinct Levy walk for each atom
    1085           0 :             idim = 0
    1086           0 :             DO ia = 1, pint_env%ndim/3
    1087           0 :                x0(1) = pint_env%x(1, 3*(ia - 1) + 1)
    1088           0 :                x0(2) = pint_env%x(1, 3*(ia - 1) + 2)
    1089           0 :                x0(3) = pint_env%x(1, 3*(ia - 1) + 3)
    1090           0 :                var = SQRT(1.0_dp/(pint_env%kT*tcorr*pint_env%mass(3*ia)))
    1091           0 :                CALL pint_levy_walk(x0, pint_env%p, var, bx, rng_gaussian)
    1092           0 :                DO ic = 1, 3
    1093           0 :                   idim = idim + 1
    1094           0 :                   DO ib = 1, pint_env%p
    1095           0 :                      pint_env%x(ib, idim) = pint_env%x(ib, idim) + bx(3*(ib - 1) + ic)
    1096             :                   END DO
    1097             :                END DO
    1098             :             END DO
    1099             : 
    1100             :          END IF
    1101             : 
    1102           0 :          DEALLOCATE (bx)
    1103           0 :          done_levy = .TRUE.
    1104             :       END IF
    1105             : 
    1106          54 :       done_init = .FALSE.
    1107          54 :       NULLIFY (input_section)
    1108             :       input_section => section_vals_get_subs_vals(pint_env%input, &
    1109          54 :                                                   "MOTION%PINT%BEADS%COORD")
    1110          54 :       CALL section_vals_get(input_section, explicit=explicit)
    1111          54 :       IF (explicit) THEN
    1112             :          CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1113           8 :                                    n_rep_val=n_rep_val)
    1114           8 :          IF (n_rep_val > 0) THEN
    1115           8 :             CPASSERT(n_rep_val == 1)
    1116             :             CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1117           8 :                                       r_vals=r_vals)
    1118           8 :             IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim) &
    1119           0 :                CPABORT("Invalid size of MOTION%PINT%BEADS%COORD")
    1120           8 :             ic = 0
    1121        9278 :             DO idim = 1, pint_env%ndim
    1122       46358 :                DO ib = 1, pint_env%p
    1123       37080 :                   ic = ic + 1
    1124       46350 :                   pint_env%x(ib, idim) = r_vals(ic)
    1125             :                END DO
    1126             :             END DO
    1127             :             done_init = .TRUE.
    1128             :          END IF
    1129             :       END IF
    1130             : 
    1131          54 :       done_rand = .FALSE.
    1132             :       CALL section_vals_val_get(pint_env%input, &
    1133             :                                 "MOTION%PINT%INIT%RANDOMIZE_POS", &
    1134          54 :                                 l_val=ltmp)
    1135          54 :       IF (ltmp) THEN
    1136             : 
    1137           0 :          IF (pint_env%beadwise_constraints) THEN
    1138             :             WRITE (UNIT=msg, FMT=*) "Beadwise constraints are not supported if "// &
    1139             :                "a random noise is applied to the initialization of the bead positions. "// &
    1140           0 :                "Please use hot start (default)."
    1141           0 :             CPABORT(msg)
    1142             :          END IF
    1143             : 
    1144           0 :          DO idim = 1, pint_env%ndim
    1145           0 :             DO ib = 1, pint_env%p
    1146             :                pint_env%x(ib, idim) = pint_env%x(ib, idim) + &
    1147             :                                       pint_env%randomG%next(variance=pint_env%beta/ &
    1148           0 :                                                             SQRT(12.0_dp*pint_env%mass(idim)))
    1149             :             END DO
    1150             :          END DO
    1151             :          done_rand = .TRUE.
    1152             :       END IF
    1153             : 
    1154          54 :       WRITE (tmp, '(A)') "Bead positions initialization:"
    1155          54 :       IF (done_init) THEN
    1156           8 :          WRITE (msg, '(A,A)') TRIM(tmp), " input structure"
    1157          46 :       ELSE IF (done_levy) THEN
    1158           0 :          WRITE (msg, '(A,A)') TRIM(tmp), " Levy random walk"
    1159             :       ELSE
    1160          46 :          WRITE (msg, '(A,A)') TRIM(tmp), " hot start"
    1161             :       END IF
    1162          54 :       CALL pint_write_line(msg)
    1163             : 
    1164          54 :       IF (done_levy) THEN
    1165           0 :          WRITE (msg, '(A,F6.3)') "Levy walk at effective temperature: ", tcorr
    1166             :       END IF
    1167             : 
    1168          54 :       IF (done_rand) THEN
    1169           0 :          WRITE (msg, '(A)') "Added gaussian noise to the positions of the beads."
    1170           0 :          CALL pint_write_line(msg)
    1171             :       END IF
    1172             : 
    1173         108 :    END SUBROUTINE pint_init_x
    1174             : 
    1175             : ! ***************************************************************************
    1176             : !> \brief  Initialize velocities
    1177             : !> \param  pint_env the pint env in which you should initialize the
    1178             : !>         velocity
    1179             : !> \par    History
    1180             : !>         2010-12-16 gathered all velocity-init code here [lwalewski]
    1181             : !>         2011-04-05 added centroid velocity initialization [lwalewski]
    1182             : !>         2011-12-19 removed optional parameter kT, target temperature is
    1183             : !>                    now determined from the input directly [lwalewski]
    1184             : !> \author fawzi
    1185             : !> \note   Initialization is done according to the following protocol:
    1186             : !>         1. set all the velocities to FORCE_EVAL%SUBSYS%VELOCITY if present
    1187             : !>         2. scale the velocities according to the actual temperature
    1188             : !>            (has no effect if vels not present in 1.)
    1189             : !>         3. draw vels for the remaining dof from MB distribution
    1190             : !>            (all or non-centroid modes only depending on 1.)
    1191             : !>         4. add random noise to the centroid vels if CENTROID_SPEED == T
    1192             : !>         5. set the vels for all dof to 0.0 if VELOCITY_QUENCH == T
    1193             : !>         6. set the vels according to the explicit values from the input
    1194             : !>            if present
    1195             : ! **************************************************************************************************
    1196          54 :    SUBROUTINE pint_init_v(pint_env)
    1197             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1198             : 
    1199             :       CHARACTER(len=default_string_length)               :: msg, stmp, stmp1, stmp2, unit_str
    1200             :       INTEGER                                            :: first_mode, i, ia, ib, ic, idim, ierr, &
    1201             :                                                             itmp, j, n_rep_val, nparticle, &
    1202             :                                                             nparticle_kind
    1203             :       LOGICAL                                            :: done_init, done_quench, done_scale, &
    1204             :                                                             done_sped, explicit, ltmp, vels_present
    1205             :       REAL(kind=dp)                                      :: actual_t, ek, factor, rtmp, target_t, &
    1206             :                                                             unit_conv
    1207          54 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: vel
    1208          54 :       REAL(kind=dp), DIMENSION(:), POINTER               :: r_vals
    1209             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1210          54 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1211             :       TYPE(cell_type), POINTER                           :: cell
    1212             :       TYPE(cp_logger_type), POINTER                      :: logger
    1213             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1214             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    1215             :       TYPE(f_env_type), POINTER                          :: f_env
    1216             :       TYPE(global_constraint_type), POINTER              :: gci
    1217             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1218          54 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    1219             :       TYPE(molecule_list_type), POINTER                  :: molecules
    1220          54 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    1221             :       TYPE(particle_list_type), POINTER                  :: particles
    1222          54 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1223             :       TYPE(section_vals_type), POINTER                   :: input_section
    1224             : 
    1225          54 :       NULLIFY (logger)
    1226         108 :       logger => cp_get_default_logger()
    1227             : 
    1228             :       ! Get constraint info, if needed
    1229             :       ! Create a force environment which will be identical to
    1230             :       ! the bead that is being processed by the processor.
    1231          54 :       IF (pint_env%simpar%constraint) THEN
    1232           6 :          NULLIFY (subsys, cell)
    1233           6 :          NULLIFY (atomic_kinds, local_particles, particles)
    1234           6 :          NULLIFY (local_molecules, molecules, molecule_kinds, gci)
    1235           6 :          NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set)
    1236             : 
    1237           6 :          CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env)
    1238           6 :          CALL force_env_get(force_env=f_env%force_env, subsys=subsys)
    1239           6 :          CALL f_env_rm_defaults(f_env, ierr)
    1240           6 :          CPASSERT(ierr == 0)
    1241             : 
    1242             :          ! Get gci and more from subsys
    1243             :          CALL cp_subsys_get(subsys=subsys, &
    1244             :                             cell=cell, &
    1245             :                             atomic_kinds=atomic_kinds, &
    1246             :                             local_particles=local_particles, &
    1247             :                             particles=particles, &
    1248             :                             local_molecules=local_molecules, &
    1249             :                             molecules=molecules, &
    1250             :                             molecule_kinds=molecule_kinds, &
    1251           6 :                             gci=gci)
    1252             : 
    1253           6 :          nparticle_kind = atomic_kinds%n_els
    1254           6 :          atomic_kind_set => atomic_kinds%els
    1255           6 :          molecule_kind_set => molecule_kinds%els
    1256           6 :          nparticle = particles%n_els
    1257           6 :          particle_set => particles%els
    1258           6 :          molecule_set => molecules%els
    1259             : 
    1260             :          ! Allocate work storage
    1261          18 :          ALLOCATE (vel(3, nparticle))
    1262          78 :          vel(:, :) = 0.0_dp
    1263             :          CALL getold(gci, local_molecules, molecule_set, &
    1264          12 :                      molecule_kind_set, particle_set, cell)
    1265             :       END IF
    1266             : 
    1267             :       ! read the velocities from the input file if they are given explicitly
    1268          54 :       vels_present = .FALSE.
    1269          54 :       NULLIFY (input_section)
    1270             :       input_section => section_vals_get_subs_vals(pint_env%input, &
    1271          54 :                                                   "FORCE_EVAL%SUBSYS%VELOCITY")
    1272          54 :       CALL section_vals_get(input_section, explicit=explicit)
    1273          54 :       IF (explicit) THEN
    1274             : 
    1275             :          CALL section_vals_val_get(input_section, "PINT_UNIT", &
    1276           2 :                                    c_val=unit_str)
    1277           2 :          unit_conv = cp_unit_to_cp2k(1.0_dp, TRIM(unit_str))
    1278             : 
    1279             :          ! assign all the beads with the same velocities from FORCE_EVAL%SUBSYS%VELOCITY
    1280           2 :          NULLIFY (r_vals)
    1281             :          CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1282           2 :                                    n_rep_val=n_rep_val)
    1283           2 :          stmp = ""
    1284           2 :          WRITE (stmp, *) n_rep_val
    1285             :          msg = "Invalid number of atoms in FORCE_EVAL%SUBSYS%VELOCITY ("// &
    1286           2 :                TRIM(ADJUSTL(stmp))//")."
    1287           2 :          IF (3*n_rep_val /= pint_env%ndim) &
    1288           0 :             CPABORT(msg)
    1289          14 :          DO ia = 1, pint_env%ndim/3
    1290             :             CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1291          12 :                                       i_rep_val=ia, r_vals=r_vals)
    1292          12 :             itmp = SIZE(r_vals)
    1293          12 :             stmp = ""
    1294          12 :             WRITE (stmp, *) itmp
    1295             :             msg = "Number of coordinates != 3 in FORCE_EVAL%SUBSYS%VELOCITY ("// &
    1296          12 :                   TRIM(ADJUSTL(stmp))//")."
    1297          12 :             IF (itmp /= 3) &
    1298           0 :                CPABORT(msg)
    1299         110 :             DO ib = 1, pint_env%p
    1300         396 :                DO ic = 1, 3
    1301         288 :                   idim = 3*(ia - 1) + ic
    1302         384 :                   pint_env%v(ib, idim) = r_vals(ic)*unit_conv
    1303             :                END DO
    1304             :             END DO
    1305             :          END DO
    1306             : 
    1307             :          vels_present = .TRUE.
    1308             :       END IF
    1309             : 
    1310             :       ! set the actual temperature...
    1311             :       IF (vels_present) THEN
    1312             :          ! ...from the initial velocities
    1313           2 :          ek = 0.0_dp
    1314          14 :          DO ia = 1, pint_env%ndim/3
    1315             :             rtmp = 0.0_dp
    1316          48 :             DO ic = 1, 3
    1317          36 :                idim = 3*(ia - 1) + ic
    1318          48 :                rtmp = rtmp + pint_env%v(1, idim)*pint_env%v(1, idim)
    1319             :             END DO
    1320          14 :             ek = ek + 0.5_dp*pint_env%mass(idim)*rtmp
    1321             :          END DO
    1322           2 :          actual_t = 2.0_dp*ek/pint_env%ndim
    1323             :       ELSE
    1324             :          ! ...using the temperature value from the input
    1325          52 :          actual_t = pint_env%kT
    1326             :       END IF
    1327             : 
    1328             :       ! set the target temperature
    1329          54 :       target_t = pint_env%kT
    1330             :       CALL section_vals_val_get(pint_env%input, &
    1331             :                                 "MOTION%PINT%INIT%VELOCITY_SCALE", &
    1332          54 :                                 l_val=done_scale)
    1333          54 :       IF (vels_present) THEN
    1334           2 :          IF (done_scale) THEN
    1335             :             ! rescale the velocities to match the target temperature
    1336           2 :             rtmp = SQRT(target_t/actual_t)
    1337          14 :             DO ia = 1, pint_env%ndim/3
    1338         110 :                DO ib = 1, pint_env%p
    1339         396 :                   DO ic = 1, 3
    1340         288 :                      idim = 3*(ia - 1) + ic
    1341         384 :                      pint_env%v(ib, idim) = rtmp*pint_env%v(ib, idim)
    1342             :                   END DO
    1343             :                END DO
    1344             :             END DO
    1345             :          ELSE
    1346             :             target_t = actual_t
    1347             :          END IF
    1348             :       END IF
    1349             : 
    1350             :       ! draw velocities from the M-B distribution...
    1351          54 :       IF (vels_present) THEN
    1352             :          ! ...for non-centroid modes only
    1353           2 :          CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
    1354           2 :          first_mode = 2
    1355             :       ELSE
    1356             :          ! ...for all the modes
    1357             :          first_mode = 1
    1358             :       END IF
    1359       64944 :       DO idim = 1, SIZE(pint_env%uv, 2)
    1360      362196 :          DO ib = first_mode, SIZE(pint_env%uv, 1)
    1361             :             pint_env%uv(ib, idim) = &
    1362      362142 :                pint_env%randomG%next(variance=target_t/pint_env%mass_fict(ib, idim))
    1363             :          END DO
    1364             :       END DO
    1365             : 
    1366             :       ! add random component to the centroid velocity if requested
    1367          54 :       done_sped = .FALSE.
    1368             :       CALL section_vals_val_get(pint_env%input, &
    1369             :                                 "MOTION%PINT%INIT%CENTROID_SPEED", &
    1370          54 :                                 l_val=ltmp)
    1371          54 :       IF (ltmp) THEN
    1372           0 :          CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
    1373           0 :          DO idim = 1, pint_env%ndim
    1374             :             rtmp = pint_env%randomG%next(variance=pint_env%mass(idim)*pint_env%kT) &
    1375           0 :                    /pint_env%mass(idim)
    1376           0 :             DO ib = 1, pint_env%p
    1377           0 :                pint_env%v(ib, idim) = pint_env%v(ib, idim) + rtmp
    1378             :             END DO
    1379             :          END DO
    1380           0 :          CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
    1381           0 :          done_sped = .TRUE.
    1382             :       END IF
    1383             : 
    1384             :       ! quench (set to zero) velocities for all the modes if requested
    1385             :       ! (disregard all the initialization done so far)
    1386          54 :       done_quench = .FALSE.
    1387             :       CALL section_vals_val_get(pint_env%input, &
    1388             :                                 "MOTION%PINT%INIT%VELOCITY_QUENCH", &
    1389          54 :                                 l_val=ltmp)
    1390          54 :       IF (ltmp) THEN
    1391           0 :          DO idim = 1, pint_env%ndim
    1392           0 :             DO ib = 1, pint_env%p
    1393           0 :                pint_env%v(ib, idim) = 0.0_dp
    1394             :             END DO
    1395             :          END DO
    1396           0 :          CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
    1397           0 :          done_quench = .TRUE.
    1398             :       END IF
    1399             : 
    1400             :       ! set the velocities to the values from the input if they are explicit
    1401             :       ! (disregard all the initialization done so far)
    1402          54 :       done_init = .FALSE.
    1403          54 :       NULLIFY (input_section)
    1404             :       input_section => section_vals_get_subs_vals(pint_env%input, &
    1405          54 :                                                   "MOTION%PINT%BEADS%VELOCITY")
    1406          54 :       CALL section_vals_get(input_section, explicit=explicit)
    1407          54 :       IF (explicit) THEN
    1408             :          CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1409           8 :                                    n_rep_val=n_rep_val)
    1410           8 :          IF (n_rep_val > 0) THEN
    1411           8 :             CPASSERT(n_rep_val == 1)
    1412             :             CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1413           8 :                                       r_vals=r_vals)
    1414           8 :             IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim) &
    1415           0 :                CPABORT("Invalid size of MOTION%PINT%BEAD%VELOCITY")
    1416           8 :             itmp = 0
    1417        9278 :             DO idim = 1, pint_env%ndim
    1418       46358 :                DO ib = 1, pint_env%p
    1419       37080 :                   itmp = itmp + 1
    1420       46350 :                   pint_env%v(ib, idim) = r_vals(itmp)
    1421             :                END DO
    1422             :             END DO
    1423           8 :             CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
    1424           8 :             done_init = .TRUE.
    1425             :          END IF
    1426             :       END IF
    1427             : 
    1428          54 :       unit_conv = cp_unit_from_cp2k(1.0_dp, "K")
    1429          54 :       WRITE (stmp1, '(F10.2)') target_t*pint_env%propagator%temp_sim2phys*unit_conv
    1430          54 :       msg = "Bead velocities initialization:"
    1431          54 :       IF (done_init) THEN
    1432           8 :          msg = TRIM(msg)//" input structure"
    1433          46 :       ELSE IF (done_quench) THEN
    1434           0 :          msg = TRIM(msg)//" quenching (set to 0.0)"
    1435             :       ELSE
    1436          46 :          IF (vels_present) THEN
    1437           2 :             msg = TRIM(ADJUSTL(msg))//" centroid +"
    1438             :          END IF
    1439          46 :          msg = TRIM(ADJUSTL(msg))//" Maxwell-Boltzmann at "//TRIM(ADJUSTL(stmp1))//" K."
    1440             :       END IF
    1441          54 :       CALL pint_write_line(msg)
    1442             : 
    1443          54 :       IF (done_init .AND. done_quench) THEN
    1444           0 :          msg = "WARNING: exclusive options requested (velocity restart and quenching)"
    1445           0 :          CPWARN(msg)
    1446           0 :          msg = "WARNING: velocity restart took precedence"
    1447           0 :          CPWARN(msg)
    1448             :       END IF
    1449             : 
    1450          54 :       IF ((.NOT. done_init) .AND. (.NOT. done_quench)) THEN
    1451          46 :          IF (vels_present .AND. done_scale) THEN
    1452           2 :             WRITE (stmp1, '(F10.2)') actual_t*unit_conv
    1453           2 :             WRITE (stmp2, '(F10.2)') target_t*unit_conv
    1454             :             msg = "Scaled initial velocities from "//TRIM(ADJUSTL(stmp1))// &
    1455           2 :                   " to "//TRIM(ADJUSTL(stmp2))//" K as requested."
    1456           2 :             CPWARN(msg)
    1457             :          END IF
    1458          46 :          IF (done_sped) THEN
    1459           0 :             msg = "Added random component to the initial centroid velocities."
    1460           0 :             CPWARN(msg)
    1461             :          END IF
    1462             :       END IF
    1463             : 
    1464             :       ! Apply constraints to the initial velocities
    1465          54 :       IF (pint_env%simpar%constraint) THEN
    1466           6 :          IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
    1467             :             ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD
    1468           0 :             factor = SQRT(REAL(pint_env%p, dp))
    1469             :          ELSE
    1470             :             ! lowest NM is centroid
    1471             :             factor = 1.0_dp
    1472             :          END IF
    1473             :          ! Beadwise constraints
    1474           6 :          IF (pint_env%beadwise_constraints) THEN
    1475           2 :             IF (pint_env%logger%para_env%is_source()) THEN
    1476           1 :                CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
    1477           5 :                DO ib = 1, pint_env%p
    1478          16 :                   DO i = 1, nparticle
    1479          52 :                      DO j = 1, 3
    1480             :                         ! Centroid is also constrained. This has to be changed if the initialization
    1481             :                         ! of the positions of the beads is done as free particles (LEVY_POS_SAMPLE)
    1482             :                         ! or if a Gaussian noise is added (RANDOMIZE_POS)
    1483          36 :                         particle_set(i)%r(j) = pint_env%x(1, j + (i - 1)*3)/factor
    1484          48 :                         vel(j, i) = pint_env%v(ib, j + (i - 1)*3)
    1485             :                      END DO
    1486             :                   END DO
    1487             :                   ! Possibly update the target values
    1488             :                   CALL shake_update_targets(gci, local_molecules, molecule_set, &
    1489             :                                             molecule_kind_set, pint_env%dt, &
    1490           4 :                                             f_env%force_env%root_section)
    1491             :                   CALL rattle_control(gci, local_molecules, molecule_set, &
    1492             :                                       molecule_kind_set, particle_set, &
    1493             :                                       vel, pint_env%dt, pint_env%simpar%shake_tol, &
    1494             :                                       pint_env%simpar%info_constraint, &
    1495             :                                       pint_env%simpar%lagrange_multipliers, &
    1496             :                                       .FALSE., &
    1497             :                                       cell, mp_comm_self, &
    1498           4 :                                       local_particles)
    1499          17 :                   DO i = 1, nparticle
    1500          52 :                      DO j = 1, 3
    1501          48 :                         pint_env%v(ib, j + (i - 1)*3) = vel(j, i)
    1502             :                      END DO
    1503             :                   END DO
    1504             :                END DO
    1505             :                ! Transform back to normal modes:
    1506           1 :                CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
    1507             :             END IF
    1508             :             ! Broadcast updated velocities to other nodes
    1509         182 :             CALL pint_env%logger%para_env%bcast(pint_env%uv)
    1510             :             ! Centroid constraints
    1511             :          ELSE
    1512             :             ! Transform positions and velocities to Cartesian coordinates:
    1513           4 :             IF (pint_env%logger%para_env%is_source()) THEN
    1514           8 :                DO i = 1, nparticle
    1515          26 :                   DO j = 1, 3
    1516          18 :                      particle_set(i)%r(j) = pint_env%x(1, j + (i - 1)*3)/factor
    1517          24 :                      vel(j, i) = pint_env%uv(1, j + (i - 1)*3)/factor
    1518             :                   END DO
    1519             :                END DO
    1520             :                ! Possibly update the target values
    1521             :                CALL shake_update_targets(gci, local_molecules, molecule_set, &
    1522             :                                          molecule_kind_set, pint_env%dt, &
    1523           2 :                                          f_env%force_env%root_section)
    1524             :                CALL rattle_control(gci, local_molecules, molecule_set, &
    1525             :                                    molecule_kind_set, particle_set, &
    1526             :                                    vel, pint_env%dt, pint_env%simpar%shake_tol, &
    1527             :                                    pint_env%simpar%info_constraint, &
    1528             :                                    pint_env%simpar%lagrange_multipliers, &
    1529             :                                    .FALSE., &
    1530             :                                    cell, mp_comm_self, &
    1531           2 :                                    local_particles)
    1532             :             END IF
    1533             :             ! Broadcast updated velocities to other nodes
    1534           4 :             CALL pint_env%logger%para_env%bcast(vel)
    1535             :             ! Transform back to normal modes
    1536          16 :             DO i = 1, nparticle
    1537          52 :                DO j = 1, 3
    1538          48 :                   pint_env%uv(1, j + (i - 1)*3) = vel(j, i)*factor
    1539             :                END DO
    1540             :             END DO
    1541             :          END IF
    1542             :       END IF
    1543             : 
    1544         108 :    END SUBROUTINE pint_init_v
    1545             : 
    1546             : ! ***************************************************************************
    1547             : !> \brief  Assign initial postions and velocities to the thermostats.
    1548             : !> \param pint_env ...
    1549             : !> \param kT ...
    1550             : !> \date   2010-12-15
    1551             : !> \author Lukasz Walewski
    1552             : !> \note   Extracted from pint_init
    1553             : ! **************************************************************************************************
    1554          54 :    SUBROUTINE pint_init_t(pint_env, kT)
    1555             : 
    1556             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1557             :       REAL(kind=dp), INTENT(in), OPTIONAL                :: kT
    1558             : 
    1559             :       INTEGER                                            :: ib, idim, ii, inos, n_rep_val
    1560             :       LOGICAL                                            :: explicit, gle_restart
    1561             :       REAL(kind=dp)                                      :: mykt
    1562          54 :       REAL(kind=dp), DIMENSION(:), POINTER               :: r_vals
    1563             :       TYPE(section_vals_type), POINTER                   :: input_section
    1564             : 
    1565         106 :       IF (pint_env%pimd_thermostat == thermostat_nose) THEN
    1566             : 
    1567          26 :          mykt = pint_env%kT
    1568          26 :          IF (PRESENT(kT)) mykt = kT
    1569        9476 :          DO idim = 1, SIZE(pint_env%tv, 3)
    1570       29096 :             DO ib = 1, SIZE(pint_env%tv, 2)
    1571       88218 :                DO inos = 1, SIZE(pint_env%tv, 1)
    1572             :                   pint_env%tv(inos, ib, idim) = &
    1573       78768 :                      pint_env%randomG%next(variance=mykt/pint_env%Q(ib))
    1574             :                END DO
    1575             :             END DO
    1576             :          END DO
    1577          26 :          IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    1578          74 :             pint_env%tv(:, 1, :) = 0.0_dp
    1579             :          END IF
    1580             : 
    1581          26 :          NULLIFY (input_section)
    1582             :          input_section => section_vals_get_subs_vals(pint_env%input, &
    1583          26 :                                                      "MOTION%PINT%NOSE%COORD")
    1584          26 :          CALL section_vals_get(input_section, explicit=explicit)
    1585          26 :          IF (explicit) THEN
    1586             :             CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1587           6 :                                       n_rep_val=n_rep_val)
    1588           6 :             IF (n_rep_val > 0) THEN
    1589           6 :                CPASSERT(n_rep_val == 1)
    1590             :                CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1591           6 :                                          r_vals=r_vals)
    1592           6 :                IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim*pint_env%nnos) &
    1593           0 :                   CPABORT("Invalid size of MOTION%PINT%NOSE%COORD")
    1594           6 :                ii = 0
    1595          60 :                DO idim = 1, pint_env%ndim
    1596         276 :                   DO ib = 1, pint_env%p
    1597         918 :                      DO inos = 1, pint_env%nnos
    1598         648 :                         ii = ii + 1
    1599         864 :                         pint_env%tx(inos, ib, idim) = r_vals(ii)
    1600             :                      END DO
    1601             :                   END DO
    1602             :                END DO
    1603             :             END IF
    1604             :          END IF
    1605          26 :          IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    1606          74 :             pint_env%tx(:, 1, :) = 0.0_dp
    1607             :          END IF
    1608             : 
    1609          26 :          NULLIFY (input_section)
    1610             :          input_section => section_vals_get_subs_vals(pint_env%input, &
    1611          26 :                                                      "MOTION%PINT%NOSE%VELOCITY")
    1612          26 :          CALL section_vals_get(input_section, explicit=explicit)
    1613          26 :          IF (explicit) THEN
    1614             :             CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1615           6 :                                       n_rep_val=n_rep_val)
    1616           6 :             IF (n_rep_val > 0) THEN
    1617           6 :                CPASSERT(n_rep_val == 1)
    1618             :                CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
    1619           6 :                                          r_vals=r_vals)
    1620           6 :                IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim*pint_env%nnos) &
    1621           0 :                   CPABORT("Invalid size of MOTION%PINT%NOSE%VELOCITY")
    1622           6 :                ii = 0
    1623          60 :                DO idim = 1, pint_env%ndim
    1624         276 :                   DO ib = 1, pint_env%p
    1625         918 :                      DO inos = 1, pint_env%nnos
    1626         648 :                         ii = ii + 1
    1627         864 :                         pint_env%tv(inos, ib, idim) = r_vals(ii)
    1628             :                      END DO
    1629             :                   END DO
    1630             :                END DO
    1631             :             END IF
    1632           6 :             IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    1633           0 :                pint_env%tv(:, 1, :) = 0.0_dp
    1634             :             END IF
    1635             :          END IF
    1636             : 
    1637          28 :       ELSEIF (pint_env%pimd_thermostat == thermostat_gle) THEN
    1638           2 :          NULLIFY (input_section)
    1639             :          input_section => section_vals_get_subs_vals(pint_env%input, &
    1640           2 :                                                      "MOTION%PINT%GLE")
    1641           2 :          CALL section_vals_get(input_section, explicit=explicit)
    1642           2 :          IF (explicit) THEN
    1643             :             CALL restart_gle(pint_env%gle, input_section, save_mem=.FALSE., &
    1644           2 :                              restart=gle_restart)
    1645             :          END IF
    1646             :       END IF
    1647             : 
    1648          54 :    END SUBROUTINE pint_init_t
    1649             : 
    1650             : ! ***************************************************************************
    1651             : !> \brief  Prepares the forces, etc. to perform an PIMD step
    1652             : !> \param pint_env ...
    1653             : !> \param helium_env ...
    1654             : !> \par    History
    1655             : !>           Added nh_energy calculation [hforbert]
    1656             : !>           Bug fixes for no thermostats [hforbert]
    1657             : !>           2016-07-14 Modified to work with independent helium_env [cschran]
    1658             : !> \author fawzi
    1659             : ! **************************************************************************************************
    1660         136 :    SUBROUTINE pint_init_f(pint_env, helium_env)
    1661             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1662             :       TYPE(helium_solvent_p_type), DIMENSION(:), &
    1663             :          OPTIONAL, POINTER                               :: helium_env
    1664             : 
    1665             :       INTEGER                                            :: ib, idim, inos
    1666             :       REAL(kind=dp)                                      :: e_h
    1667             :       TYPE(cp_logger_type), POINTER                      :: logger
    1668             : 
    1669          68 :       NULLIFY (logger)
    1670          68 :       logger => cp_get_default_logger()
    1671             : 
    1672             :       ! initialize iteration info
    1673          68 :       CALL cp_iterate(logger%iter_info, iter_nr=pint_env%first_step)
    1674          68 :       CALL cp_iterate(pint_env%logger%iter_info, iter_nr=pint_env%first_step)
    1675             : 
    1676          68 :       CALL pint_x2u(pint_env)
    1677          68 :       CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h)
    1678          68 :       CALL pint_calc_f(pint_env)
    1679             : 
    1680             :       ! add helium forces to the solute's internal ones
    1681             :       ! Assume that helium has been already initialized and helium_env(1)
    1682             :       ! contains proper forces in force_avrg array at ionode
    1683          68 :       IF (PRESENT(helium_env)) THEN
    1684          14 :          IF (logger%para_env%is_source()) THEN
    1685         574 :             pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
    1686             :          END IF
    1687        2282 :          CALL logger%para_env%bcast(pint_env%f)
    1688             :       END IF
    1689          68 :       CALL pint_f2uf(pint_env)
    1690             : 
    1691             :       ! set the centroid forces to 0 if FIX_CENTROID_POS
    1692          68 :       IF (pint_env%first_propagated_mode .EQ. 2) THEN
    1693           0 :          pint_env%uf(1, :) = 0.0_dp
    1694             :       END IF
    1695             : 
    1696          68 :       CALL pint_calc_e_kin_beads_u(pint_env)
    1697          68 :       CALL pint_calc_e_vir(pint_env)
    1698       65084 :       DO idim = 1, SIZE(pint_env%uf_h, 2)
    1699      363380 :          DO ib = pint_env%first_propagated_mode, SIZE(pint_env%uf_h, 1)
    1700      363312 :             pint_env%uf(ib, idim) = REAL(pint_env%nrespa, dp)*pint_env%uf(ib, idim)
    1701             :          END DO
    1702             :       END DO
    1703             : 
    1704          68 :       IF (pint_env%nnos > 0) THEN
    1705        9576 :          DO idim = 1, SIZE(pint_env%uf_h, 2)
    1706       29556 :             DO ib = 1, SIZE(pint_env%uf_h, 1)
    1707             :                pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* &
    1708       29520 :                                            pint_env%uv(ib, idim)**2 - pint_env%kT)/pint_env%Q(ib)
    1709             :             END DO
    1710             :          END DO
    1711             : 
    1712        9576 :          DO idim = 1, pint_env%ndim
    1713       29556 :             DO ib = 1, pint_env%p
    1714       60228 :                DO inos = 1, pint_env%nnos - 1
    1715             :                   pint_env%tf(inos + 1, ib, idim) = pint_env%tv(inos, ib, idim)**2 - &
    1716       60228 :                                                     pint_env%kT/pint_env%Q(ib)
    1717             :                END DO
    1718       69768 :                DO inos = 1, pint_env%nnos - 1
    1719             :                   pint_env%tf(inos, ib, idim) = pint_env%tf(inos, ib, idim) &
    1720       60228 :                                                 - pint_env%tv(inos, ib, idim)*pint_env%tv(inos + 1, ib, idim)
    1721             :                END DO
    1722             :             END DO
    1723             :          END DO
    1724          36 :          CALL pint_calc_nh_energy(pint_env)
    1725             :       END IF
    1726             : 
    1727          68 :    END SUBROUTINE pint_init_f
    1728             : 
    1729             : ! ***************************************************************************
    1730             : !> \brief  Perform the PIMD simulation (main MD loop)
    1731             : !> \param pint_env ...
    1732             : !> \param globenv ...
    1733             : !> \param helium_env ...
    1734             : !> \par    History
    1735             : !>         2003-11 created [fawzi]
    1736             : !>         renamed from pint_run to pint_do_run because of conflicting name
    1737             : !>           of pint_run in input_constants [hforbert]
    1738             : !>         2009-12-14 globenv parameter added to handle soft exit
    1739             : !>           requests [lwalewski]
    1740             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1741             : !> \author Fawzi Mohamed
    1742             : !> \note   Everything should be read for an md step.
    1743             : ! **************************************************************************************************
    1744          54 :    SUBROUTINE pint_do_run(pint_env, globenv, helium_env)
    1745             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1746             :       TYPE(global_environment_type), POINTER             :: globenv
    1747             :       TYPE(helium_solvent_p_type), DIMENSION(:), &
    1748             :          OPTIONAL, POINTER                               :: helium_env
    1749             : 
    1750             :       INTEGER                                            :: k, step
    1751             :       LOGICAL                                            :: should_stop
    1752             :       REAL(kind=dp)                                      :: scal
    1753             :       TYPE(cp_logger_type), POINTER                      :: logger
    1754             :       TYPE(f_env_type), POINTER                          :: f_env
    1755             : 
    1756             :       ! initialize iteration info
    1757          54 :       CALL cp_iterate(pint_env%logger%iter_info, iter_nr=pint_env%first_step)
    1758             : 
    1759             :       ! iterate replica pint counter by accessing the globally saved
    1760             :       ! force environment error/logger variables and setting them
    1761             :       ! explicitly to the pimd "PINT" step value
    1762             :       CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
    1763          54 :                               f_env=f_env)
    1764          54 :       NULLIFY (logger)
    1765          54 :       logger => cp_get_default_logger()
    1766             :       CALL cp_iterate(logger%iter_info, &
    1767          54 :                       iter_nr=pint_env%first_step)
    1768          54 :       CALL f_env_rm_defaults(f_env)
    1769             : 
    1770          54 :       pint_env%iter = pint_env%first_step
    1771             : 
    1772          54 :       IF (PRESENT(helium_env)) THEN
    1773          14 :          IF (ASSOCIATED(helium_env)) THEN
    1774             :             ! set the properties accumulated over the whole MC process to 0
    1775          34 :             DO k = 1, SIZE(helium_env)
    1776          80 :                helium_env(k)%helium%proarea%accu(:) = 0.0_dp
    1777          80 :                helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
    1778          80 :                helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
    1779          80 :                helium_env(k)%helium%mominer%accu(:) = 0.0_dp
    1780          20 :                IF (helium_env(k)%helium%rho_present) THEN
    1781    40202020 :                   helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
    1782             :                END IF
    1783          34 :                IF (helium_env(k)%helium%rdf_present) THEN
    1784       20020 :                   helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
    1785             :                END IF
    1786             :             END DO
    1787             :          END IF
    1788             :       END IF
    1789             : 
    1790             :       ! write the properties at 0-th step
    1791          54 :       CALL pint_calc_energy(pint_env)
    1792          54 :       CALL pint_calc_total_action(pint_env)
    1793          54 :       CALL pint_write_ener(pint_env)
    1794          54 :       CALL pint_write_action(pint_env)
    1795          54 :       CALL pint_write_centroids(pint_env)
    1796          54 :       CALL pint_write_trajectory(pint_env)
    1797          54 :       CALL pint_write_com(pint_env)
    1798          54 :       CALL pint_write_rgyr(pint_env)
    1799             : 
    1800             :       ! main PIMD loop
    1801         482 :       DO step = 1, pint_env%num_steps
    1802             : 
    1803         428 :          pint_env%iter = pint_env%iter + 1
    1804             :          CALL cp_iterate(pint_env%logger%iter_info, &
    1805             :                          last=(step == pint_env%num_steps), &
    1806         428 :                          iter_nr=pint_env%iter)
    1807             :          CALL cp_iterate(logger%iter_info, &
    1808             :                          last=(step == pint_env%num_steps), &
    1809         428 :                          iter_nr=pint_env%iter)
    1810         428 :          pint_env%t = pint_env%t + pint_env%dt
    1811             : 
    1812         428 :          IF (pint_env%t_tol > 0.0_dp) THEN
    1813           0 :             IF (ABS(2._dp*pint_env%e_kin_beads/(pint_env%p*pint_env%ndim) &
    1814             :                     - pint_env%kT) > pint_env%t_tol) THEN
    1815           0 :                scal = SQRT(pint_env%kT*(pint_env%p*pint_env%ndim)/(2.0_dp*pint_env%e_kin_beads))
    1816           0 :                pint_env%uv = scal*pint_env%uv
    1817           0 :                CALL pint_init_f(pint_env, helium_env=helium_env)
    1818             :             END IF
    1819             :          END IF
    1820         794 :          CALL pint_step(pint_env, helium_env=helium_env)
    1821             : 
    1822         428 :          CALL pint_write_ener(pint_env)
    1823         428 :          CALL pint_write_action(pint_env)
    1824         428 :          CALL pint_write_centroids(pint_env)
    1825         428 :          CALL pint_write_trajectory(pint_env)
    1826         428 :          CALL pint_write_com(pint_env)
    1827         428 :          CALL pint_write_rgyr(pint_env)
    1828             : 
    1829             :          CALL write_restart(root_section=pint_env%input, &
    1830         428 :                             pint_env=pint_env, helium_env=helium_env)
    1831             : 
    1832             :          ! exit from the main loop if soft exit has been requested
    1833         428 :          CALL external_control(should_stop, "PINT", globenv=globenv)
    1834         482 :          IF (should_stop) EXIT
    1835             : 
    1836             :       END DO
    1837             : 
    1838             :       ! remove iteration level
    1839          54 :       CALL cp_rm_iter_level(pint_env%logger%iter_info, "PINT")
    1840             : 
    1841          54 :    END SUBROUTINE pint_do_run
    1842             : 
    1843             : ! ***************************************************************************
    1844             : !> \brief  Performs a scan of the helium-solute interaction energy
    1845             : !> \param pint_env ...
    1846             : !> \param helium_env ...
    1847             : !> \date   2013-11-26
    1848             : !> \parm   History
    1849             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1850             : !> \author Lukasz Walewski
    1851             : ! **************************************************************************************************
    1852           0 :    SUBROUTINE pint_run_scan(pint_env, helium_env)
    1853             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1854             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1855             : 
    1856             :       CHARACTER(len=default_string_length)               :: comment
    1857             :       INTEGER                                            :: unit_nr
    1858           0 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: DATA
    1859             :       TYPE(section_vals_type), POINTER                   :: print_key
    1860             : 
    1861           0 :       NULLIFY (pint_env%logger, print_key)
    1862           0 :       pint_env%logger => cp_get_default_logger()
    1863             : 
    1864             :       ! assume that ionode always has at least one helium_env
    1865           0 :       IF (pint_env%logger%para_env%is_source()) THEN
    1866             :          print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1867           0 :                                                  "MOTION%PINT%HELIUM%PRINT%RHO")
    1868             :       END IF
    1869             : 
    1870             :       ! perform the actual scan wrt the COM of the solute
    1871           0 :       CALL helium_intpot_scan(pint_env, helium_env)
    1872             : 
    1873             :       ! output the interaction potential into a cubefile
    1874             :       ! assume that ionode always has at least one helium_env
    1875           0 :       IF (pint_env%logger%para_env%is_source()) THEN
    1876             : 
    1877             :          unit_nr = cp_print_key_unit_nr( &
    1878             :                    pint_env%logger, &
    1879             :                    print_key, &
    1880             :                    middle_name="helium-pot", &
    1881             :                    extension=".cube", &
    1882             :                    file_position="REWIND", &
    1883           0 :                    do_backup=.FALSE.)
    1884             : 
    1885           0 :          comment = "Solute - helium interaction potential"
    1886           0 :          NULLIFY (DATA)
    1887           0 :          DATA => helium_env(1)%helium%rho_inst(1, :, :, :)
    1888             :          CALL helium_write_cubefile( &
    1889             :             unit_nr, &
    1890             :             comment, &
    1891             :             helium_env(1)%helium%center - 0.5_dp* &
    1892             :             (helium_env(1)%helium%rho_maxr - helium_env(1)%helium%rho_delr), &
    1893             :             helium_env(1)%helium%rho_delr, &
    1894             :             helium_env(1)%helium%rho_nbin, &
    1895           0 :             DATA)
    1896             : 
    1897           0 :          CALL m_flush(unit_nr)
    1898           0 :          CALL cp_print_key_finished_output(unit_nr, pint_env%logger, print_key)
    1899             : 
    1900             :       END IF
    1901             : 
    1902             :       ! output solute positions
    1903           0 :       CALL pint_write_centroids(pint_env)
    1904           0 :       CALL pint_write_trajectory(pint_env)
    1905             : 
    1906           0 :    END SUBROUTINE pint_run_scan
    1907             : 
    1908             : ! ***************************************************************************
    1909             : !> \brief  Does an PINT step (and nrespa harmonic evaluations)
    1910             : !> \param pint_env ...
    1911             : !> \param helium_env ...
    1912             : !> \par    History
    1913             : !>           various bug fixes [hforbert]
    1914             : !>           10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl]
    1915             : !>           04.2016 Changed to work with helium_env [cschran]
    1916             : !>           10.2018 Added centroid constraints [cschran+rperez]
    1917             : !>           10.2021 Added beadwise constraints [lduran]
    1918             : !> \author fawzi
    1919             : ! **************************************************************************************************
    1920         428 :    SUBROUTINE pint_step(pint_env, helium_env)
    1921             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    1922             :       TYPE(helium_solvent_p_type), DIMENSION(:), &
    1923             :          OPTIONAL, POINTER                               :: helium_env
    1924             : 
    1925             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pint_step'
    1926             : 
    1927             :       INTEGER                                            :: handle, i, ia, ib, idim, ierr, inos, &
    1928             :                                                             iresp, j, k, nbeads, nparticle, &
    1929             :                                                             nparticle_kind
    1930             :       REAL(kind=dp)                                      :: dt_temp, dti, dti2, dti22, e_h, factor, &
    1931             :                                                             rn, tdti, time_start, time_stop, tol
    1932         428 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: pos, vel
    1933         428 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: tmp
    1934             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1935         428 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1936             :       TYPE(cell_type), POINTER                           :: cell
    1937             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1938             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    1939             :       TYPE(f_env_type), POINTER                          :: f_env
    1940             :       TYPE(global_constraint_type), POINTER              :: gci
    1941             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1942         428 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    1943             :       TYPE(molecule_list_type), POINTER                  :: molecules
    1944         428 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    1945             :       TYPE(particle_list_type), POINTER                  :: particles
    1946         428 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1947             : 
    1948         428 :       CALL timeset(routineN, handle)
    1949         428 :       time_start = m_walltime()
    1950             : 
    1951         428 :       rn = REAL(pint_env%nrespa, dp)
    1952         428 :       dti = pint_env%dt/rn
    1953         428 :       dti2 = dti/2._dp
    1954         428 :       tdti = 2.*dti
    1955         428 :       dti22 = dti**2/2._dp
    1956             : 
    1957             :       ! Get constraint info, if needed
    1958             :       ! Create a force environment which will be identical to
    1959             :       ! the bead that is being processed by the processor.
    1960         428 :       IF (pint_env%simpar%constraint) THEN
    1961          24 :          NULLIFY (subsys, cell)
    1962          24 :          NULLIFY (atomic_kinds, local_particles, particles)
    1963          24 :          NULLIFY (local_molecules, molecules, molecule_kinds, gci)
    1964          24 :          NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set)
    1965             : 
    1966          24 :          CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env)
    1967          24 :          CALL force_env_get(force_env=f_env%force_env, subsys=subsys)
    1968          24 :          CALL f_env_rm_defaults(f_env, ierr)
    1969          24 :          CPASSERT(ierr == 0)
    1970             : 
    1971             :          ! Get gci and more from subsys
    1972             :          CALL cp_subsys_get(subsys=subsys, &
    1973             :                             cell=cell, &
    1974             :                             atomic_kinds=atomic_kinds, &
    1975             :                             local_particles=local_particles, &
    1976             :                             particles=particles, &
    1977             :                             local_molecules=local_molecules, &
    1978             :                             molecules=molecules, &
    1979             :                             molecule_kinds=molecule_kinds, &
    1980          24 :                             gci=gci)
    1981             : 
    1982          24 :          nparticle_kind = atomic_kinds%n_els
    1983          24 :          atomic_kind_set => atomic_kinds%els
    1984          24 :          molecule_kind_set => molecule_kinds%els
    1985          24 :          nparticle = particles%n_els
    1986          24 :          nbeads = pint_env%p
    1987          24 :          particle_set => particles%els
    1988          24 :          molecule_set => molecules%els
    1989             : 
    1990             :          ! Allocate work storage
    1991          72 :          ALLOCATE (pos(3, nparticle))
    1992          72 :          ALLOCATE (vel(3, nparticle))
    1993         312 :          pos(:, :) = 0.0_dp
    1994         312 :          vel(:, :) = 0.0_dp
    1995             : 
    1996          24 :          IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
    1997             :             ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD
    1998           0 :             factor = SQRT(REAL(pint_env%p, dp))
    1999             :          ELSE
    2000             :             factor = 1.0_dp
    2001             :          END IF
    2002             : 
    2003             :          CALL getold(gci, local_molecules, molecule_set, &
    2004          48 :                      molecule_kind_set, particle_set, cell)
    2005             :       END IF
    2006             : 
    2007         690 :       SELECT CASE (pint_env%harm_integrator)
    2008             :       CASE (integrate_numeric)
    2009             : 
    2010         938 :          DO iresp = 1, pint_env%nrespa
    2011             : 
    2012             :             ! integrate bead positions, first_propagated_mode = { 1, 2 }
    2013             :             ! Nose needs an extra step
    2014         676 :             IF (pint_env%pimd_thermostat == thermostat_nose) THEN
    2015             : 
    2016             :                !Set thermostat action of constrained DoF to zero:
    2017         616 :                IF (pint_env%simpar%constraint) THEN
    2018          48 :                   DO k = 1, pint_env%n_atoms_constraints
    2019          32 :                      ia = pint_env%atoms_constraints(k)
    2020         144 :                      DO j = 3*(ia - 1) + 1, 3*ia
    2021         416 :                         pint_env%tv(:, 1, j) = 0.0_dp
    2022             :                      END DO
    2023             :                   END DO
    2024             :                END IF
    2025             : 
    2026             :                ! Exempt centroid from thermostat for CMD
    2027         616 :                IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2028        5920 :                   pint_env%tx(:, 1, :) = 0.0_dp
    2029        5920 :                   pint_env%tv(:, 1, :) = 0.0_dp
    2030        5920 :                   pint_env%tf(:, 1, :) = 0.0_dp
    2031             :                END IF
    2032             : 
    2033        4832 :                DO i = pint_env%first_propagated_mode, pint_env%p
    2034             :                   pint_env%ux(i, :) = pint_env%ux(i, :) - &
    2035       93968 :                                       dti22*pint_env%uv(i, :)*pint_env%tv(1, i, :)
    2036             :                END DO
    2037      411700 :                pint_env%tx = pint_env%tx + dti*pint_env%tv + dti22*pint_env%tf
    2038             : 
    2039         616 :                IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2040        5920 :                   pint_env%tx(:, 1, :) = 0.0_dp
    2041        5920 :                   pint_env%tv(:, 1, :) = 0.0_dp
    2042        5920 :                   pint_env%tf(:, 1, :) = 0.0_dp
    2043             :                END IF
    2044             : 
    2045             :             END IF
    2046             :             !Integrate position in harmonic springs (uf_h) and physical potential
    2047             :             !(uf)
    2048        5124 :             DO i = pint_env%first_propagated_mode, pint_env%p
    2049             :                pint_env%ux_t(i, :) = pint_env%ux(i, :) + &
    2050             :                                      dti*pint_env%uv(i, :) + &
    2051             :                                      dti22*(pint_env%uf_h(i, :) + &
    2052      133140 :                                             pint_env%uf(i, :))
    2053             :             END DO
    2054             : 
    2055             :             ! apply thermostats to velocities
    2056        1292 :             SELECT CASE (pint_env%pimd_thermostat)
    2057             :             CASE (thermostat_nose)
    2058             : 
    2059         616 :                IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2060        5920 :                   pint_env%tx(:, 1, :) = 0.0_dp
    2061        5920 :                   pint_env%tv(:, 1, :) = 0.0_dp
    2062        5920 :                   pint_env%tf(:, 1, :) = 0.0_dp
    2063             :                END IF
    2064             : 
    2065             :                pint_env%uv_t = pint_env%uv - dti2* &
    2066      230368 :                                pint_env%uv*pint_env%tv(1, :, :)
    2067         616 :                tmp => pint_env%tv_t
    2068         616 :                pint_env%tv_t => pint_env%tv
    2069         616 :                pint_env%tv => tmp
    2070      411700 :                pint_env%tv = pint_env%tv_old + tdti*pint_env%tf
    2071      411700 :                pint_env%tv_old = pint_env%tv_t
    2072      411700 :                pint_env%tv_t = pint_env%tv_t + dti2*pint_env%tf
    2073             :             CASE DEFAULT
    2074       58492 :                pint_env%uv_t = pint_env%uv
    2075             :             END SELECT
    2076             : 
    2077             :             !Set thermostat action of constrained DoF to zero:
    2078         676 :             IF (pint_env%simpar%constraint) THEN
    2079          48 :                DO k = 1, pint_env%n_atoms_constraints
    2080          32 :                   ia = pint_env%atoms_constraints(k)
    2081         144 :                   DO j = 3*(ia - 1) + 1, 3*ia
    2082         384 :                      pint_env%tv(:, 1, j) = 0.0_dp
    2083         416 :                      pint_env%tv_t(:, 1, j) = 0.0_dp
    2084             :                   END DO
    2085             :                END DO
    2086             :             END IF
    2087             : 
    2088             :             ! Exempt centroid from thermostat for CMD
    2089         676 :             IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2090        5920 :                pint_env%tx(:, 1, :) = 0.0_dp
    2091        5920 :                pint_env%tv(:, 1, :) = 0.0_dp
    2092        5920 :                pint_env%tf(:, 1, :) = 0.0_dp
    2093             :             END IF
    2094             : 
    2095             :             !Integrate harmonic velocities and physical velocities
    2096      173368 :             pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf)
    2097             : 
    2098             :             ! physical forces are only applied in first respa step.
    2099      173368 :             pint_env%uf = 0.0_dp
    2100             :             ! calc harmonic forces at new pos
    2101      173368 :             pint_env%ux = pint_env%ux_t
    2102             : 
    2103             :             ! Apply centroid constraints (SHAKE)
    2104         676 :             IF (pint_env%simpar%constraint) THEN
    2105          16 :                IF (pint_env%logger%para_env%is_source()) THEN
    2106          32 :                   DO i = 1, nparticle
    2107         104 :                      DO j = 1, 3
    2108          72 :                         pos(j, i) = pint_env%ux(1, j + (i - 1)*3)
    2109          96 :                         vel(j, i) = pint_env%uv_t(1, j + (i - 1)*3)
    2110             :                      END DO
    2111             :                   END DO
    2112             : 
    2113             :                   ! Possibly update the target values
    2114             :                   CALL shake_update_targets(gci, local_molecules, molecule_set, &
    2115             :                                             molecule_kind_set, dti, &
    2116           8 :                                             f_env%force_env%root_section)
    2117             :                   CALL shake_control(gci, local_molecules, molecule_set, &
    2118             :                                      molecule_kind_set, particle_set, &
    2119             :                                      pos, vel, dti, pint_env%simpar%shake_tol, &
    2120             :                                      pint_env%simpar%info_constraint, &
    2121             :                                      pint_env%simpar%lagrange_multipliers, &
    2122             :                                      pint_env%simpar%dump_lm, cell, &
    2123           8 :                                      mp_comm_self, local_particles)
    2124             :                END IF
    2125             :                ! Positions and velocities of centroid were constrained by SHAKE
    2126          16 :                CALL pint_env%logger%para_env%bcast(pos)
    2127          16 :                CALL pint_env%logger%para_env%bcast(vel)
    2128             :                ! Transform back to normal modes:
    2129          64 :                DO i = 1, nparticle
    2130         208 :                   DO j = 1, 3
    2131         144 :                      pint_env%ux(1, j + (i - 1)*3) = pos(j, i)
    2132         192 :                      pint_env%uv_t(1, j + (i - 1)*3) = vel(j, i)
    2133             :                   END DO
    2134             :                END DO
    2135             : 
    2136             :             END IF
    2137             :             ! Exempt centroid from thermostat for CMD
    2138         676 :             IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2139        5920 :                pint_env%tx(:, 1, :) = 0.0_dp
    2140        5920 :                pint_env%tv(:, 1, :) = 0.0_dp
    2141        5920 :                pint_env%tf(:, 1, :) = 0.0_dp
    2142             :             END IF
    2143             : 
    2144         676 :             CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h)
    2145      173368 :             pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf)
    2146             : 
    2147             :             ! For last respa step include integration of physical and helium
    2148             :             ! forces
    2149         676 :             IF (iresp == pint_env%nrespa) THEN
    2150         262 :                CALL pint_u2x(pint_env)
    2151         262 :                CALL pint_calc_f(pint_env)
    2152             :                ! perform helium step and add helium forces
    2153         262 :                IF (PRESENT(helium_env)) THEN
    2154          50 :                   CALL helium_step(helium_env, pint_env)
    2155             :                   !Update force of solute in pint_env
    2156          50 :                   IF (pint_env%logger%para_env%is_source()) THEN
    2157        1150 :                      pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
    2158             :                   END IF
    2159        4550 :                   CALL pint_env%logger%para_env%bcast(pint_env%f)
    2160             :                END IF
    2161             : 
    2162         262 :                CALL pint_f2uf(pint_env)
    2163             :                ! set the centroid forces to 0 if FIX_CENTROID_POS
    2164         262 :                IF (pint_env%first_propagated_mode .EQ. 2) THEN
    2165           0 :                   pint_env%uf(1, :) = 0.0_dp
    2166             :                END IF
    2167             :                !Scale physical forces and integrate velocities with physical
    2168             :                !forces
    2169      128944 :                pint_env%uf = pint_env%uf*rn
    2170      128944 :                pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf
    2171             : 
    2172             :             END IF
    2173             : 
    2174             :             ! Apply second half of thermostats
    2175         938 :             SELECT CASE (pint_env%pimd_thermostat)
    2176             :             CASE (thermostat_nose)
    2177             :                ! Exempt centroid from thermostat for CMD
    2178         616 :                IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2179        5920 :                   pint_env%tx(:, 1, :) = 0.0_dp
    2180        5920 :                   pint_env%tv(:, 1, :) = 0.0_dp
    2181        5920 :                   pint_env%tf(:, 1, :) = 0.0_dp
    2182             :                END IF
    2183        3758 :                DO i = 1, 6
    2184        3380 :                   tol = 0._dp
    2185     1353452 :                   pint_env%uv_new = pint_env%uv_t/(1.+dti2*pint_env%tv(1, :, :))
    2186      154976 :                   DO idim = 1, pint_env%ndim
    2187      678416 :                      DO ib = 1, pint_env%p
    2188             :                         pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* &
    2189             :                                                     pint_env%uv_new(ib, idim)**2 - pint_env%kT*pint_env%kTcorr)/ &
    2190      675036 :                                                    pint_env%Q(ib)
    2191             :                      END DO
    2192             :                   END DO
    2193             : 
    2194             :                   !Set thermostat action of constrained DoF to zero:
    2195        3380 :                   IF (pint_env%simpar%constraint) THEN
    2196         288 :                      DO k = 1, pint_env%n_atoms_constraints
    2197         192 :                         ia = pint_env%atoms_constraints(k)
    2198         864 :                         DO j = 3*(ia - 1) + 1, 3*ia
    2199        2496 :                            pint_env%tf(:, 1, j) = 0.0_dp
    2200             :                         END DO
    2201             :                      END DO
    2202             :                   END IF
    2203             : 
    2204             :                   ! Exempt centroid from thermostat for CMD
    2205        3380 :                   IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2206       35520 :                      pint_env%tx(:, 1, :) = 0.0_dp
    2207       35520 :                      pint_env%tv(:, 1, :) = 0.0_dp
    2208       35520 :                      pint_env%tf(:, 1, :) = 0.0_dp
    2209             :                   END IF
    2210             : 
    2211      154976 :                   DO idim = 1, pint_env%ndim
    2212      678416 :                      DO ib = 1, pint_env%p
    2213     1743120 :                         DO inos = 1, pint_env%nnos - 1
    2214             :                            pint_env%tv_new(inos, ib, idim) = &
    2215             :                               (pint_env%tv_t(inos, ib, idim) + dti2*pint_env%tf(inos, ib, idim))/ &
    2216     1219680 :                               (1._dp + dti2*pint_env%tv(inos + 1, ib, idim))
    2217             :                            pint_env%tf(inos + 1, ib, idim) = &
    2218             :                               (pint_env%tv_new(inos, ib, idim)**2 - &
    2219     1219680 :                                pint_env%kT*pint_env%kTcorr/pint_env%Q(ib))
    2220             :                            tol = MAX(tol, ABS(pint_env%tv(inos, ib, idim) &
    2221     1743120 :                                               - pint_env%tv_new(inos, ib, idim)))
    2222             :                         END DO
    2223             :                         !Set thermostat action of constrained DoF to zero:
    2224      523440 :                         IF (pint_env%simpar%constraint) THEN
    2225       10368 :                            DO k = 1, pint_env%n_atoms_constraints
    2226        6912 :                               ia = pint_env%atoms_constraints(k)
    2227       31104 :                               DO j = 3*(ia - 1) + 1, 3*ia
    2228       82944 :                                  pint_env%tv_new(:, 1, j) = 0.0_dp
    2229       89856 :                                  pint_env%tf(:, 1, j) = 0.0_dp
    2230             :                               END DO
    2231             :                            END DO
    2232             :                         END IF
    2233             : 
    2234             :                         ! Exempt centroid from thermostat for CMD
    2235      523440 :                         IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2236     3196800 :                            pint_env%tx(:, 1, :) = 0.0_dp
    2237     3196800 :                            pint_env%tv(:, 1, :) = 0.0_dp
    2238     3196800 :                            pint_env%tf(:, 1, :) = 0.0_dp
    2239             :                         END IF
    2240             : 
    2241             :                         pint_env%tv_new(pint_env%nnos, ib, idim) = &
    2242             :                            pint_env%tv_t(pint_env%nnos, ib, idim) + &
    2243      523440 :                            dti2*pint_env%tf(pint_env%nnos, ib, idim)
    2244             :                         tol = MAX(tol, ABS(pint_env%tv(pint_env%nnos, ib, idim) &
    2245      523440 :                                            - pint_env%tv_new(pint_env%nnos, ib, idim)))
    2246             :                         tol = MAX(tol, ABS(pint_env%uv(ib, idim) &
    2247      523440 :                                            - pint_env%uv_new(ib, idim)))
    2248             :                         !Set thermostat action of constrained DoF to zero:
    2249      523440 :                         IF (pint_env%simpar%constraint) THEN
    2250       10368 :                            DO k = 1, pint_env%n_atoms_constraints
    2251        6912 :                               ia = pint_env%atoms_constraints(k)
    2252       31104 :                               DO j = 3*(ia - 1) + 1, 3*ia
    2253       89856 :                                  pint_env%tv_new(:, 1, j) = 0.0_dp
    2254             :                               END DO
    2255             :                            END DO
    2256             :                         END IF
    2257             :                         ! Exempt centroid from thermostat for CMD
    2258      675036 :                         IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2259     3196800 :                            pint_env%tx(:, 1, :) = 0.0_dp
    2260     3196800 :                            pint_env%tv(:, 1, :) = 0.0_dp
    2261     3196800 :                            pint_env%tf(:, 1, :) = 0.0_dp
    2262             :                         END IF
    2263             : 
    2264             :                      END DO
    2265             :                   END DO
    2266             : 
    2267      678416 :                   pint_env%uv = pint_env%uv_new
    2268     2421536 :                   pint_env%tv = pint_env%tv_new
    2269        3380 :                   IF (tol <= pint_env%v_tol) EXIT
    2270             :                   ! Exempt centroid from thermostat for CMD
    2271        3758 :                   IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2272       35520 :                      pint_env%tx(:, 1, :) = 0.0_dp
    2273       35520 :                      pint_env%tv(:, 1, :) = 0.0_dp
    2274       35520 :                      pint_env%tf(:, 1, :) = 0.0_dp
    2275             :                   END IF
    2276             :                END DO
    2277             : 
    2278             :                ! Apply centroid constraints (RATTLE)
    2279         616 :                IF (pint_env%simpar%constraint) THEN
    2280          16 :                   IF (pint_env%logger%para_env%is_source()) THEN
    2281             :                      ! Reset particle r, due to force calc:
    2282          32 :                      DO i = 1, nparticle
    2283         104 :                         DO j = 1, 3
    2284          72 :                            vel(j, i) = pint_env%uv(1, j + (i - 1)*3)
    2285          96 :                            particle_set(i)%r(j) = pint_env%ux(1, j + (i - 1)*3)
    2286             :                         END DO
    2287             :                      END DO
    2288             : 
    2289             :                      ! Small time step for all small integrations steps
    2290             :                      ! Big step for last RESPA
    2291           8 :                      IF (iresp == pint_env%nrespa) THEN
    2292           4 :                         dt_temp = dti
    2293             :                      ELSE
    2294           4 :                         dt_temp = dti*rn
    2295             :                      END IF
    2296             :                      CALL rattle_control(gci, local_molecules, molecule_set, &
    2297             :                                          molecule_kind_set, particle_set, &
    2298             :                                          vel, dt_temp, pint_env%simpar%shake_tol, &
    2299             :                                          pint_env%simpar%info_constraint, &
    2300             :                                          pint_env%simpar%lagrange_multipliers, &
    2301             :                                          pint_env%simpar%dump_lm, cell, &
    2302           8 :                                          mp_comm_self, local_particles)
    2303             :                   END IF
    2304             :                   ! Velocities of centroid were constrained by RATTLE
    2305             :                   ! Broadcast updated velocities to other nodes
    2306          16 :                   CALL pint_env%logger%para_env%bcast(vel)
    2307             : 
    2308          64 :                   DO i = 1, nparticle
    2309         208 :                      DO j = 1, 3
    2310         192 :                         pint_env%uv(1, j + (i - 1)*3) = vel(j, i)
    2311             :                      END DO
    2312             :                   END DO
    2313             :                END IF
    2314             : 
    2315        2048 :                DO inos = 1, pint_env%nnos - 1
    2316             :                   pint_env%tf(inos, :, :) = pint_env%tf(inos, :, :) &
    2317      264200 :                                             - pint_env%tv(inos, :, :)*pint_env%tv(inos + 1, :, :)
    2318             :                END DO
    2319             : 
    2320             :                ! Exempt centroid from thermostat for CMD
    2321         616 :                IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
    2322        5920 :                   pint_env%tx(:, 1, :) = 0.0_dp
    2323        5920 :                   pint_env%tv(:, 1, :) = 0.0_dp
    2324        5920 :                   pint_env%tf(:, 1, :) = 0.0_dp
    2325             :                END IF
    2326             : 
    2327             :             CASE (thermostat_gle)
    2328           4 :                CALL pint_gle_step(pint_env)
    2329       55300 :                pint_env%uv = pint_env%uv_t
    2330             :             CASE DEFAULT
    2331        3196 :                pint_env%uv = pint_env%uv_t
    2332             :             END SELECT
    2333             :          END DO
    2334             : 
    2335             :       CASE (integrate_exact)
    2336             :          ! The Liouvillian splitting is as follows:
    2337             :          ! 1. Thermostat
    2338             :          ! 2. 0.5*physical integration
    2339             :          ! 3. Exact harmonic integration + apply constraints (SHAKE)
    2340             :          ! 4. 0.5*physical integration
    2341             :          ! 5. Thermostat + apply constraints (RATTLE)
    2342             : 
    2343             :          ! 1. Apply thermostats
    2344         268 :          SELECT CASE (pint_env%pimd_thermostat)
    2345             :          CASE (thermostat_pile)
    2346             :             CALL pint_pile_step(vold=pint_env%uv, &
    2347             :                                 vnew=pint_env%uv_t, &
    2348             :                                 p=pint_env%p, &
    2349             :                                 ndim=pint_env%ndim, &
    2350             :                                 first_mode=pint_env%first_propagated_mode, &
    2351             :                                 masses=pint_env%mass_fict, &
    2352         102 :                                 pile_therm=pint_env%pile_therm)
    2353             :          CASE (thermostat_piglet)
    2354             :             CALL pint_piglet_step(vold=pint_env%uv, &
    2355             :                                   vnew=pint_env%uv_t, &
    2356             :                                   first_mode=pint_env%first_propagated_mode, &
    2357             :                                   masses=pint_env%mass_fict, &
    2358          10 :                                   piglet_therm=pint_env%piglet_therm)
    2359             :          CASE (thermostat_qtb)
    2360             :             CALL pint_qtb_step(vold=pint_env%uv, &
    2361             :                                vnew=pint_env%uv_t, &
    2362             :                                p=pint_env%p, &
    2363             :                                ndim=pint_env%ndim, &
    2364             :                                masses=pint_env%mass_fict, &
    2365          30 :                                qtb_therm=pint_env%qtb_therm)
    2366             :          CASE DEFAULT
    2367        1246 :             pint_env%uv_t = pint_env%uv
    2368             :          END SELECT
    2369             : 
    2370             :          ! 2. 1/2*Physical integration
    2371     1531858 :          pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf
    2372             : 
    2373             :          ! 3. Exact harmonic integration
    2374         166 :          IF (pint_env%first_propagated_mode == 1) THEN
    2375             :             ! The centroid is integrated via standard velocity-verlet
    2376             :             ! Commented out code is only there to show similarities to
    2377             :             ! Numeric integrator
    2378             :             pint_env%ux_t(1, :) = pint_env%ux(1, :) + &
    2379      231610 :                                   dti*pint_env%uv_t(1, :) !+ &
    2380             :             !                      dti22*pint_env%uf_h(1, :)
    2381             :             !pint_env%uv_t(1, :) = pint_env%uv_t(1, :)+ &
    2382             :             !                      dti2*pint_env%uf_h(1, :)
    2383             :          ELSE
    2384             :             ! set velocities to zero for fixed centroids
    2385           0 :             pint_env%ux_t(1, :) = pint_env%ux(1, :)
    2386           0 :             pint_env%uv_t(1, :) = 0.0_dp
    2387             :          END IF
    2388             :          ! Other modes are integrated exactly
    2389        1392 :          DO i = 2, pint_env%p
    2390             :             pint_env%ux_t(i, :) = pint_env%cosex(i)*pint_env%ux(i, :) &
    2391     1070030 :                                   + pint_env%iwsinex(i)*pint_env%uv_t(i, :)
    2392             :             pint_env%uv_t(i, :) = pint_env%cosex(i)*pint_env%uv_t(i, :) &
    2393     1070196 :                                   - pint_env%wsinex(i)*pint_env%ux(i, :)
    2394             :          END DO
    2395             : 
    2396             :          ! Apply constraints (SHAKE)
    2397         166 :          IF (pint_env%simpar%constraint) THEN
    2398             :             ! Beadwise constraints
    2399          16 :             IF (pint_env%beadwise_constraints) THEN
    2400           8 :                IF (pint_env%logger%para_env%is_source()) THEN
    2401             :                   ! Transform positions and velocities to Cartesian coordinates:
    2402           4 :                   CALL pint_u2x(pint_env, ux=pint_env%ux_t, x=pint_env%x)
    2403           4 :                   CALL pint_u2x(pint_env, ux=pint_env%uv_t, x=pint_env%v)
    2404          20 :                   DO ib = 1, nbeads
    2405          64 :                      DO i = 1, nparticle
    2406         208 :                         DO j = 1, 3
    2407         144 :                            pos(j, i) = pint_env%x(ib, j + (i - 1)*3)
    2408         192 :                            vel(j, i) = pint_env%v(ib, j + (i - 1)*3)
    2409             :                         END DO
    2410             :                      END DO
    2411             :                      ! Possibly update the target values
    2412             :                      CALL shake_update_targets(gci, local_molecules, molecule_set, &
    2413             :                                                molecule_kind_set, dti, &
    2414          16 :                                                f_env%force_env%root_section)
    2415             :                      CALL shake_control(gci, local_molecules, molecule_set, &
    2416             :                                         molecule_kind_set, particle_set, &
    2417             :                                         pos, vel, dti, pint_env%simpar%shake_tol, &
    2418             :                                         pint_env%simpar%info_constraint, &
    2419             :                                         pint_env%simpar%lagrange_multipliers, &
    2420             :                                         pint_env%simpar%dump_lm, cell, &
    2421          16 :                                         mp_comm_self, local_particles)
    2422          68 :                      DO i = 1, nparticle
    2423         208 :                         DO j = 1, 3
    2424         144 :                            pint_env%x(ib, j + (i - 1)*3) = pos(j, i)
    2425         192 :                            pint_env%v(ib, j + (i - 1)*3) = vel(j, i)
    2426             :                         END DO
    2427             :                      END DO
    2428             :                   END DO
    2429             :                   ! Transform back to normal modes:
    2430           4 :                   CALL pint_x2u(pint_env, ux=pint_env%ux_t, x=pint_env%x)
    2431           4 :                   CALL pint_x2u(pint_env, ux=pint_env%uv_t, x=pint_env%v)
    2432             :                END IF
    2433             :                ! Broadcast positions and velocities to all nodes
    2434         728 :                CALL pint_env%logger%para_env%bcast(pint_env%ux_t)
    2435         728 :                CALL pint_env%logger%para_env%bcast(pint_env%uv_t)
    2436             :                ! Centroid constraints
    2437             :             ELSE
    2438           8 :                IF (pint_env%logger%para_env%is_source()) THEN
    2439             :                   ! Transform positions and velocities to Cartesian coordinates:
    2440          16 :                   DO i = 1, nparticle
    2441          52 :                      DO j = 1, 3
    2442          36 :                         pos(j, i) = pint_env%ux_t(1, j + (i - 1)*3)/factor
    2443          48 :                         vel(j, i) = pint_env%uv_t(1, j + (i - 1)*3)/factor
    2444             :                      END DO
    2445             :                   END DO
    2446             :                   ! Possibly update the target values
    2447             :                   CALL shake_update_targets(gci, local_molecules, molecule_set, &
    2448             :                                             molecule_kind_set, dti, &
    2449           4 :                                             f_env%force_env%root_section)
    2450             :                   CALL shake_control(gci, local_molecules, molecule_set, &
    2451             :                                      molecule_kind_set, particle_set, &
    2452             :                                      pos, vel, dti, pint_env%simpar%shake_tol, &
    2453             :                                      pint_env%simpar%info_constraint, &
    2454             :                                      pint_env%simpar%lagrange_multipliers, &
    2455             :                                      pint_env%simpar%dump_lm, cell, &
    2456           4 :                                      mp_comm_self, local_particles)
    2457             :                END IF
    2458             :                ! Broadcast positions and velocities to all nodes
    2459           8 :                CALL pint_env%logger%para_env%bcast(pos)
    2460           8 :                CALL pint_env%logger%para_env%bcast(vel)
    2461             :                ! Transform back to normal modes:
    2462          32 :                DO i = 1, nparticle
    2463         104 :                   DO j = 1, 3
    2464          72 :                      pint_env%ux_t(1, j + (i - 1)*3) = pos(j, i)*factor
    2465          96 :                      pint_env%uv_t(1, j + (i - 1)*3) = vel(j, i)*factor
    2466             :                   END DO
    2467             :                END DO
    2468             :             END IF
    2469             :             ! Positions and velocities were constrained by SHAKE
    2470             :          END IF
    2471             :          ! Update positions
    2472     1531858 :          pint_env%ux = pint_env%ux_t
    2473             : 
    2474             :          ! 4. 1/2*Physical integration
    2475     1531858 :          pint_env%uf = 0.0_dp
    2476         166 :          CALL pint_u2x(pint_env)
    2477         166 :          CALL pint_calc_f(pint_env)
    2478             :          ! perform helium step and add helium forces
    2479         166 :          IF (PRESENT(helium_env)) THEN
    2480          12 :             CALL helium_step(helium_env, pint_env)
    2481             :             !Update force of solute in pint_env
    2482          12 :             IF (pint_env%logger%para_env%is_source()) THEN
    2483        1032 :                pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
    2484             :             END IF
    2485        4116 :             CALL pint_env%logger%para_env%bcast(pint_env%f)
    2486             :          END IF
    2487         166 :          CALL pint_f2uf(pint_env)
    2488             :          ! set the centroid forces to 0 if FIX_CENTROID_POS
    2489         166 :          IF (pint_env%first_propagated_mode .EQ. 2) THEN
    2490           0 :             pint_env%uf(1, :) = 0.0_dp
    2491             :          END IF
    2492     1531858 :          pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf
    2493             : 
    2494             :          ! 5. Apply thermostats
    2495         268 :          SELECT CASE (pint_env%pimd_thermostat)
    2496             :          CASE (thermostat_pile)
    2497             :             CALL pint_pile_step(vold=pint_env%uv_t, &
    2498             :                                 vnew=pint_env%uv, &
    2499             :                                 p=pint_env%p, &
    2500             :                                 ndim=pint_env%ndim, &
    2501             :                                 first_mode=pint_env%first_propagated_mode, &
    2502             :                                 masses=pint_env%mass_fict, &
    2503         102 :                                 pile_therm=pint_env%pile_therm)
    2504             :          CASE (thermostat_piglet)
    2505             :             CALL pint_piglet_step(vold=pint_env%uv_t, &
    2506             :                                   vnew=pint_env%uv, &
    2507             :                                   first_mode=pint_env%first_propagated_mode, &
    2508             :                                   masses=pint_env%mass_fict, &
    2509          10 :                                   piglet_therm=pint_env%piglet_therm)
    2510             :          CASE (thermostat_qtb)
    2511             :             CALL pint_qtb_step(vold=pint_env%uv_t, &
    2512             :                                vnew=pint_env%uv, &
    2513             :                                p=pint_env%p, &
    2514             :                                ndim=pint_env%ndim, &
    2515             :                                masses=pint_env%mass_fict, &
    2516          30 :                                qtb_therm=pint_env%qtb_therm)
    2517             :          CASE DEFAULT
    2518        1246 :             pint_env%uv = pint_env%uv_t
    2519             :          END SELECT
    2520             : 
    2521             :          ! Apply constraints (RATTLE)
    2522         594 :          IF (pint_env%simpar%constraint) THEN
    2523             :             ! Beadwise constraints
    2524          16 :             IF (pint_env%beadwise_constraints) THEN
    2525           8 :                IF (pint_env%logger%para_env%is_source()) THEN
    2526             :                   ! Transform positions and velocities to Cartesian coordinates:
    2527             :                   ! Reset particle r, due to force calc:
    2528           4 :                   CALL pint_u2x(pint_env, ux=pint_env%ux, x=pint_env%x)
    2529           4 :                   CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
    2530          20 :                   DO ib = 1, nbeads
    2531          64 :                      DO i = 1, nparticle
    2532         208 :                         DO j = 1, 3
    2533         144 :                            particle_set(i)%r(j) = pint_env%x(ib, j + (i - 1)*3)
    2534         192 :                            vel(j, i) = pint_env%v(ib, j + (i - 1)*3)
    2535             :                         END DO
    2536             :                      END DO
    2537             :                      CALL rattle_control(gci, local_molecules, &
    2538             :                                          molecule_set, molecule_kind_set, &
    2539             :                                          particle_set, vel, dti, &
    2540             :                                          pint_env%simpar%shake_tol, &
    2541             :                                          pint_env%simpar%info_constraint, &
    2542             :                                          pint_env%simpar%lagrange_multipliers, &
    2543             :                                          pint_env%simpar%dump_lm, cell, &
    2544          16 :                                          mp_comm_self, local_particles)
    2545          68 :                      DO i = 1, nparticle
    2546         208 :                         DO j = 1, 3
    2547         192 :                            pint_env%v(ib, j + (i - 1)*3) = vel(j, i)
    2548             :                         END DO
    2549             :                      END DO
    2550             :                   END DO
    2551             :                   ! Transform back to normal modes:
    2552           4 :                   CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
    2553             :                END IF
    2554         728 :                CALL pint_env%logger%para_env%bcast(pint_env%uv)
    2555             :                ! Centroid constraints
    2556             :             ELSE
    2557           8 :                IF (pint_env%logger%para_env%is_source()) THEN
    2558             :                   ! Transform positions and velocities to Cartesian coordinates:
    2559             :                   ! Reset particle r, due to force calc:
    2560          16 :                   DO i = 1, nparticle
    2561          52 :                      DO j = 1, 3
    2562          36 :                         vel(j, i) = pint_env%uv(1, j + (i - 1)*3)/factor
    2563          48 :                         particle_set(i)%r(j) = pint_env%ux(1, j + (i - 1)*3)/factor
    2564             :                      END DO
    2565             :                   END DO
    2566             :                   CALL rattle_control(gci, local_molecules, &
    2567             :                                       molecule_set, molecule_kind_set, &
    2568             :                                       particle_set, vel, dti, &
    2569             :                                       pint_env%simpar%shake_tol, &
    2570             :                                       pint_env%simpar%info_constraint, &
    2571             :                                       pint_env%simpar%lagrange_multipliers, &
    2572             :                                       pint_env%simpar%dump_lm, cell, &
    2573           4 :                                       mp_comm_self, local_particles)
    2574             :                END IF
    2575             :                ! Velocities of centroid were constrained by RATTLE
    2576             :                ! Broadcast updated velocities to other nodes
    2577           8 :                CALL pint_env%logger%para_env%bcast(vel)
    2578             : 
    2579             :                ! Transform back to normal modes:
    2580             :                ! Multiply with SQRT(n_beads) due to normal mode transformation
    2581          32 :                DO i = 1, nparticle
    2582         104 :                   DO j = 1, 3
    2583          96 :                      pint_env%uv(1, j + (i - 1)*3) = vel(j, i)*factor
    2584             :                   END DO
    2585             :                END DO
    2586             :             END IF
    2587             :          END IF
    2588             : 
    2589             :       END SELECT
    2590             : 
    2591         428 :       IF (pint_env%simpar%constraint) THEN
    2592          24 :          DEALLOCATE (pos, vel)
    2593             :       END IF
    2594             : 
    2595             :       ! calculate the energy components
    2596         428 :       CALL pint_calc_energy(pint_env)
    2597         428 :       CALL pint_calc_total_action(pint_env)
    2598             : 
    2599             :       ! check that the number of PINT steps matches
    2600             :       ! the number of force evaluations done so far
    2601             : !TODO make this check valid if we start from ITERATION != 0
    2602             : !     CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id,&
    2603             : !          f_env=f_env,new_error=new_error)
    2604             : !     NULLIFY(logger)
    2605             : !     logger => cp_error_get_logger(new_error)
    2606             : !     IF(logger%iter_info%iteration(2)/=pint_env%iter+1)&
    2607             : !        CPABORT("md & force_eval lost sychro")
    2608             : !     CALL f_env_rm_defaults(f_env,new_error,ierr)
    2609             : 
    2610         428 :       time_stop = m_walltime()
    2611         428 :       pint_env%time_per_step = time_stop - time_start
    2612         428 :       CALL pint_write_step_info(pint_env)
    2613         428 :       CALL timestop(handle)
    2614             : 
    2615         856 :    END SUBROUTINE pint_step
    2616             : 
    2617             : ! ***************************************************************************
    2618             : !> \brief  Calculate the energy components (private wrapper function)
    2619             : !> \param pint_env ...
    2620             : !> \date   2011-01-07
    2621             : !> \author Lukasz Walewski
    2622             : ! **************************************************************************************************
    2623         964 :    SUBROUTINE pint_calc_energy(pint_env)
    2624             : 
    2625             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    2626             : 
    2627             :       REAL(KIND=dp)                                      :: e_h
    2628             : 
    2629         482 :       CALL pint_calc_e_kin_beads_u(pint_env)
    2630         482 :       CALL pint_calc_e_vir(pint_env)
    2631             : 
    2632         482 :       CALL pint_calc_uf_h(pint_env, e_h=e_h)
    2633         482 :       pint_env%e_pot_h = e_h
    2634             : 
    2635         738 :       SELECT CASE (pint_env%pimd_thermostat)
    2636             :       CASE (thermostat_nose)
    2637         256 :          CALL pint_calc_nh_energy(pint_env)
    2638             :       CASE (thermostat_gle)
    2639           6 :          CALL pint_calc_gle_energy(pint_env)
    2640             :       CASE (thermostat_pile)
    2641         110 :          CALL pint_calc_pile_energy(pint_env)
    2642             :       CASE (thermostat_qtb)
    2643          36 :          CALL pint_calc_qtb_energy(pint_env)
    2644             :       CASE (thermostat_piglet)
    2645         482 :          CALL pint_calc_piglet_energy(pint_env)
    2646             :       END SELECT
    2647             : 
    2648             :       pint_env%energy(e_kin_thermo_id) = &
    2649             :          (0.5_dp*REAL(pint_env%p, dp)*REAL(pint_env%ndim, dp)*pint_env%kT - &
    2650         482 :           pint_env%e_pot_h)*pint_env%propagator%temp_sim2phys
    2651             : 
    2652        3778 :       pint_env%energy(e_potential_id) = SUM(pint_env%e_pot_bead)
    2653             : 
    2654             :       pint_env%energy(e_conserved_id) = &
    2655             :          pint_env%energy(e_potential_id)*pint_env%propagator%physpotscale + &
    2656             :          pint_env%e_pot_h + &
    2657             :          pint_env%e_kin_beads + &
    2658             :          pint_env%e_pot_t + &
    2659             :          pint_env%e_kin_t + &
    2660         482 :          pint_env%e_gle + pint_env%e_pile + pint_env%e_piglet + pint_env%e_qtb
    2661             : 
    2662             :       pint_env%energy(e_potential_id) = &
    2663         482 :          pint_env%energy(e_potential_id)/REAL(pint_env%p, dp)
    2664             : 
    2665         482 :    END SUBROUTINE pint_calc_energy
    2666             : 
    2667             : ! ***************************************************************************
    2668             : !> \brief  Calculate the harmonic force in the u basis
    2669             : !> \param  pint_env the path integral environment in which the harmonic
    2670             : !>         forces should be calculated
    2671             : !> \param e_h ...
    2672             : !> \par    History
    2673             : !>           Added normal mode transformation [hforbert]
    2674             : !> \author fawzi
    2675             : ! **************************************************************************************************
    2676        1226 :    SUBROUTINE pint_calc_uf_h(pint_env, e_h)
    2677             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    2678             :       REAL(KIND=dp), INTENT(OUT)                         :: e_h
    2679             : 
    2680        1226 :       IF (pint_env%transform == transformation_stage) THEN
    2681             :          CALL staging_calc_uf_h(pint_env%staging_env, &
    2682             :                                 pint_env%mass_beads, &
    2683             :                                 pint_env%ux, &
    2684             :                                 pint_env%uf_h, &
    2685           0 :                                 pint_env%e_pot_h)
    2686             :       ELSE
    2687             :          CALL normalmode_calc_uf_h(pint_env%normalmode_env, &
    2688             :                                    pint_env%mass_beads, &
    2689             :                                    pint_env%ux, &
    2690             :                                    pint_env%uf_h, &
    2691        1226 :                                    pint_env%e_pot_h)
    2692             :       END IF
    2693        1226 :       e_h = pint_env%e_pot_h
    2694     2559782 :       pint_env%uf_h = pint_env%uf_h/pint_env%mass_fict
    2695        1226 :    END SUBROUTINE pint_calc_uf_h
    2696             : 
    2697             : ! ***************************************************************************
    2698             : !> \brief calculates the force (and energy) in each bead, returns the sum
    2699             : !>      of the potential energy
    2700             : !> \param pint_env path integral environment on which you want to calculate
    2701             : !>        the forces
    2702             : !> \param x positions at which you want to evaluate the forces
    2703             : !> \param f the forces
    2704             : !> \param e potential energy on each bead
    2705             : !> \par    History
    2706             : !>           2009-06-15 moved helium calls out from here [lwalewski]
    2707             : !> \author fawzi
    2708             : ! **************************************************************************************************
    2709         496 :    SUBROUTINE pint_calc_f(pint_env, x, f, e)
    2710             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2711             :       REAL(kind=dp), DIMENSION(:, :), INTENT(in), &
    2712             :          OPTIONAL, TARGET                                :: x
    2713             :       REAL(kind=dp), DIMENSION(:, :), INTENT(out), &
    2714             :          OPTIONAL, TARGET                                :: f
    2715             :       REAL(kind=dp), DIMENSION(:), INTENT(out), &
    2716             :          OPTIONAL, TARGET                                :: e
    2717             : 
    2718             :       INTEGER                                            :: ib, idim
    2719         496 :       REAL(kind=dp), DIMENSION(:), POINTER               :: my_e
    2720         496 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: my_f, my_x
    2721             : 
    2722         496 :       my_x => pint_env%x
    2723           0 :       IF (PRESENT(x)) my_x => x
    2724         496 :       my_f => pint_env%f
    2725         496 :       IF (PRESENT(f)) my_f => f
    2726         496 :       my_e => pint_env%e_pot_bead
    2727         496 :       IF (PRESENT(e)) my_e => e
    2728      336286 :       DO idim = 1, pint_env%ndim
    2729     2024182 :          DO ib = 1, pint_env%p
    2730     2023686 :             pint_env%replicas%r(idim, ib) = my_x(ib, idim)
    2731             :          END DO
    2732             :       END DO
    2733         496 :       CALL rep_env_calc_e_f(pint_env%replicas, calc_f=.TRUE.)
    2734      336286 :       DO idim = 1, pint_env%ndim
    2735     2024182 :          DO ib = 1, pint_env%p
    2736             :             !ljw: is that fine ? - idim <-> ib
    2737     2023686 :             my_f(ib, idim) = pint_env%replicas%f(idim, ib)
    2738             :          END DO
    2739             :       END DO
    2740        3904 :       my_e = pint_env%replicas%f(SIZE(pint_env%replicas%f, 1), :)
    2741             : 
    2742         496 :    END SUBROUTINE pint_calc_f
    2743             : 
    2744             : ! ***************************************************************************
    2745             : !> \brief  Calculate the kinetic energy of the beads (in the u variables)
    2746             : !> \param pint_env ...
    2747             : !> \param uv ...
    2748             : !> \param e_k ...
    2749             : !> \par    History
    2750             : !>         Bug fix to give my_uv a default location if not given in call [hforbert]
    2751             : !> \author fawzi
    2752             : ! **************************************************************************************************
    2753         550 :    SUBROUTINE pint_calc_e_kin_beads_u(pint_env, uv, e_k)
    2754             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    2755             :       REAL(kind=dp), DIMENSION(:, :), INTENT(in), &
    2756             :          OPTIONAL, TARGET                                :: uv
    2757             :       REAL(kind=dp), INTENT(out), OPTIONAL               :: e_k
    2758             : 
    2759             :       INTEGER                                            :: ib, idim
    2760             :       REAL(kind=dp)                                      :: res
    2761         550 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: my_uv
    2762             : 
    2763         550 :       res = -1.0_dp
    2764         550 :       my_uv => pint_env%uv
    2765           0 :       IF (PRESENT(uv)) my_uv => uv
    2766         550 :       res = 0._dp
    2767      401230 :       DO idim = 1, pint_env%ndim
    2768     2386414 :          DO ib = 1, pint_env%p
    2769     2385864 :             res = res + pint_env%mass_fict(ib, idim)*my_uv(ib, idim)**2
    2770             :          END DO
    2771             :       END DO
    2772         550 :       res = res*0.5
    2773         550 :       IF (.NOT. PRESENT(uv)) pint_env%e_kin_beads = res
    2774         550 :       IF (PRESENT(e_k)) e_k = res
    2775         550 :    END SUBROUTINE pint_calc_e_kin_beads_u
    2776             : 
    2777             : ! ***************************************************************************
    2778             : !> \brief  Calculate the virial estimator of the real (quantum) kinetic energy
    2779             : !> \param pint_env ...
    2780             : !> \param e_vir ...
    2781             : !> \author hforbert
    2782             : !> \note   This subroutine modifies pint_env%energy(e_kin_virial_id) global
    2783             : !>         variable [lwalewski]
    2784             : ! **************************************************************************************************
    2785         550 :    ELEMENTAL SUBROUTINE pint_calc_e_vir(pint_env, e_vir)
    2786             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    2787             :       REAL(kind=dp), INTENT(out), OPTIONAL               :: e_vir
    2788             : 
    2789             :       INTEGER                                            :: ib, idim
    2790             :       REAL(kind=dp)                                      :: res, xcentroid
    2791             : 
    2792             :       res = -1.0_dp
    2793         550 :       res = 0._dp
    2794      401230 :       DO idim = 1, pint_env%ndim
    2795             :          ! calculate the centroid
    2796      400680 :          xcentroid = 0._dp
    2797     2385864 :          DO ib = 1, pint_env%p
    2798     2385864 :             xcentroid = xcentroid + pint_env%x(ib, idim)
    2799             :          END DO
    2800      400680 :          xcentroid = xcentroid/REAL(pint_env%p, dp)
    2801     2386414 :          DO ib = 1, pint_env%p
    2802     2385864 :             res = res + (pint_env%x(ib, idim) - xcentroid)*pint_env%f(ib, idim)
    2803             :          END DO
    2804             :       END DO
    2805             :       res = 0.5_dp*(REAL(pint_env%ndim, dp)* &
    2806         550 :                     (pint_env%kT*pint_env%propagator%temp_sim2phys) - res/REAL(pint_env%p, dp))
    2807         550 :       pint_env%energy(e_kin_virial_id) = res
    2808         550 :       IF (PRESENT(e_vir)) e_vir = res
    2809         550 :    END SUBROUTINE pint_calc_e_vir
    2810             : 
    2811             : ! ***************************************************************************
    2812             : !> \brief calculates the energy (potential and kinetic) of the Nose-Hoover
    2813             : !>      chain thermostats
    2814             : !> \param pint_env the path integral environment
    2815             : !> \author fawzi
    2816             : ! **************************************************************************************************
    2817         292 :    ELEMENTAL SUBROUTINE pint_calc_nh_energy(pint_env)
    2818             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    2819             : 
    2820             :       INTEGER                                            :: ib, idim, inos
    2821             :       REAL(kind=dp)                                      :: ekin, epot
    2822             : 
    2823         292 :       ekin = 0._dp
    2824       39928 :       DO idim = 1, pint_env%ndim
    2825      131008 :          DO ib = 1, pint_env%p
    2826      407412 :             DO inos = 1, pint_env%nnos
    2827      367776 :                ekin = ekin + pint_env%Q(ib)*pint_env%tv(inos, ib, idim)**2
    2828             :             END DO
    2829             :          END DO
    2830             :       END DO
    2831         292 :       pint_env%e_kin_t = 0.5_dp*ekin
    2832         292 :       epot = 0._dp
    2833       39928 :       DO idim = 1, pint_env%ndim
    2834      131008 :          DO ib = 1, pint_env%p
    2835      407412 :             DO inos = 1, pint_env%nnos
    2836      367776 :                epot = epot + pint_env%tx(inos, ib, idim)
    2837             :             END DO
    2838             :          END DO
    2839             :       END DO
    2840         292 :       pint_env%e_pot_t = pint_env%kT*epot
    2841         292 :    END SUBROUTINE pint_calc_nh_energy
    2842             : 
    2843             : ! ***************************************************************************
    2844             : !> \brief calculates the total link action of the PI system (excluding helium)
    2845             : !> \param pint_env the path integral environment
    2846             : !> \return ...
    2847             : !> \author Felix Uhl
    2848             : ! **************************************************************************************************
    2849         482 :    ELEMENTAL FUNCTION pint_calc_total_link_action(pint_env) RESULT(link_action)
    2850             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2851             :       REAL(KIND=dp)                                      :: link_action
    2852             : 
    2853             :       INTEGER                                            :: iatom, ibead, idim, indx
    2854             :       REAL(KIND=dp)                                      :: hb2m, tau, tmp_link_action
    2855             :       REAL(KIND=dp), DIMENSION(3)                        :: r
    2856             : 
    2857             :       !tau = 1/(k_B T p)
    2858         482 :       tau = pint_env%beta/REAL(pint_env%p, dp)
    2859             : 
    2860         482 :       link_action = 0.0_dp
    2861      112370 :       DO iatom = 1, pint_env%ndim/3
    2862             :          ! hbar / (2.0*m)
    2863      111888 :          hb2m = 1.0_dp/pint_env%mass((iatom - 1)*3 + 1)
    2864      111888 :          tmp_link_action = 0.0_dp
    2865      562296 :          DO ibead = 1, pint_env%p - 1
    2866     1801632 :             DO idim = 1, 3
    2867     1351224 :                indx = (iatom - 1)*3 + idim
    2868     1801632 :                r(idim) = pint_env%x(ibead, indx) - pint_env%x(ibead + 1, indx)
    2869             :             END DO
    2870      562296 :             tmp_link_action = tmp_link_action + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
    2871             :          END DO
    2872      447552 :          DO idim = 1, 3
    2873      335664 :             indx = (iatom - 1)*3 + idim
    2874      447552 :             r(idim) = pint_env%x(pint_env%p, indx) - pint_env%x(1, indx)
    2875             :          END DO
    2876      111888 :          tmp_link_action = tmp_link_action + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
    2877      112370 :          link_action = link_action + tmp_link_action/hb2m
    2878             :       END DO
    2879             : 
    2880         482 :       link_action = link_action/(2.0_dp*tau)
    2881             : 
    2882         482 :    END FUNCTION pint_calc_total_link_action
    2883             : 
    2884             : ! ***************************************************************************
    2885             : !> \brief calculates the potential action of the PI system (excluding helium)
    2886             : !> \param pint_env the path integral environment
    2887             : !> \return ...
    2888             : !> \author Felix Uhl
    2889             : ! **************************************************************************************************
    2890         482 :    ELEMENTAL FUNCTION pint_calc_total_pot_action(pint_env) RESULT(pot_action)
    2891             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2892             :       REAL(KIND=dp)                                      :: pot_action
    2893             : 
    2894             :       REAL(KIND=dp)                                      :: tau
    2895             : 
    2896         482 :       tau = pint_env%beta/REAL(pint_env%p, dp)
    2897        3778 :       pot_action = tau*SUM(pint_env%e_pot_bead)
    2898             : 
    2899         482 :    END FUNCTION pint_calc_total_pot_action
    2900             : 
    2901             : ! ***************************************************************************
    2902             : !> \brief calculates the total action of the PI system (excluding helium)
    2903             : !> \param pint_env the path integral environment
    2904             : !> \author Felix Uhl
    2905             : ! **************************************************************************************************
    2906         482 :    ELEMENTAL SUBROUTINE pint_calc_total_action(pint_env)
    2907             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
    2908             : 
    2909         482 :       pint_env%pot_action = pint_calc_total_pot_action(pint_env)
    2910         482 :       pint_env%link_action = pint_calc_total_link_action(pint_env)
    2911             : 
    2912         482 :    END SUBROUTINE pint_calc_total_action
    2913             : 
    2914           2 : END MODULE pint_methods

Generated by: LCOV version 1.15