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

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

Generated by: LCOV version 1.15