LCOV - code coverage report
Current view: top level - src/motion - helium_sampling.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 565 652 86.7 %
Date: 2024-04-20 06:29:22 Functions: 6 6 100.0 %

          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 for sampling helium variables
      10             : !> \author Lukasz Walewski
      11             : !> \date   2009-06-10
      12             : ! **************************************************************************************************
      13             : MODULE helium_sampling
      14             : 
      15             :    USE cp_external_control,             ONLY: external_control
      16             :    USE cp_log_handling,                 ONLY: cp_add_default_logger,&
      17             :                                               cp_get_default_logger,&
      18             :                                               cp_logger_create,&
      19             :                                               cp_logger_release,&
      20             :                                               cp_logger_type,&
      21             :                                               cp_rm_default_logger
      22             :    USE cp_output_handling,              ONLY: cp_add_iter_level,&
      23             :                                               cp_iterate,&
      24             :                                               cp_rm_iter_level
      25             :    USE global_types,                    ONLY: global_environment_type
      26             :    USE helium_common,                   ONLY: &
      27             :         helium_boxmean_3d, helium_calc_plength, helium_calc_rdf, helium_calc_rho, helium_com, &
      28             :         helium_eval_expansion, helium_pbc, helium_rotate, helium_set_rdf_coord_system, &
      29             :         helium_spline, helium_total_moment_of_inertia, helium_total_projected_area, &
      30             :         helium_total_winding_number, helium_update_transition_matrix
      31             :    USE helium_interactions,             ONLY: helium_bead_solute_e_f,&
      32             :                                               helium_calc_energy,&
      33             :                                               helium_solute_e_f
      34             :    USE helium_io,                       ONLY: &
      35             :         helium_print_accepts, helium_print_action, helium_print_coordinates, helium_print_energy, &
      36             :         helium_print_force, helium_print_perm, helium_print_plength, helium_print_rdf, &
      37             :         helium_print_rho, helium_print_vector, helium_write_line
      38             :    USE helium_types,                    ONLY: e_id_total,&
      39             :                                               helium_solvent_p_type,&
      40             :                                               helium_solvent_type
      41             :    USE helium_worm,                     ONLY: helium_sample_worm
      42             :    USE input_constants,                 ONLY: &
      43             :         helium_forces_average, helium_forces_last, helium_mdist_exponential, &
      44             :         helium_mdist_gaussian, helium_mdist_linear, helium_mdist_quadratic, helium_mdist_singlev, &
      45             :         helium_mdist_uniform, helium_sampling_ceperley, helium_sampling_worm
      46             :    USE input_cp2k_restarts,             ONLY: write_restart
      47             :    USE kinds,                           ONLY: default_string_length,&
      48             :                                               dp
      49             :    USE machine,                         ONLY: m_walltime
      50             :    USE physcon,                         ONLY: angstrom
      51             :    USE pint_public,                     ONLY: pint_com_pos
      52             :    USE pint_types,                      ONLY: pint_env_type
      53             :    USE splines_types,                   ONLY: spline_data_type
      54             : #include "../base/base_uses.f90"
      55             : 
      56             :    IMPLICIT NONE
      57             : 
      58             :    PRIVATE
      59             : 
      60             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      61             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_sampling'
      62             : 
      63             :    PUBLIC :: helium_do_run
      64             :    PUBLIC :: helium_sample
      65             :    PUBLIC :: helium_step
      66             : 
      67             : CONTAINS
      68             : 
      69             : ! ***************************************************************************
      70             : !> \brief  Performs MC simulation for helium (only)
      71             : !> \param helium_env ...
      72             : !> \param globenv ...
      73             : !> \date   2009-07-14
      74             : !> \par    History
      75             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
      76             : !> \author Lukasz Walewski
      77             : !> \note   This routine gets called only when HELIUM_ONLY is set to .TRUE.,
      78             : !>         so do not put any property calculations here!
      79             : ! **************************************************************************************************
      80          10 :    SUBROUTINE helium_do_run(helium_env, globenv)
      81             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
      82             :       TYPE(global_environment_type), POINTER             :: globenv
      83             : 
      84             :       INTEGER                                            :: k, num_steps, step, tot_steps
      85             :       LOGICAL                                            :: should_stop
      86             :       TYPE(cp_logger_type), POINTER                      :: logger
      87             :       TYPE(pint_env_type)                                :: pint_env
      88             : 
      89          10 :       NULLIFY (logger)
      90          20 :       logger => cp_get_default_logger()
      91             : 
      92          10 :       num_steps = 0
      93             : 
      94          10 :       IF (ASSOCIATED(helium_env)) THEN
      95          20 :          DO k = 1, SIZE(helium_env)
      96             : 
      97          10 :             NULLIFY (helium_env(k)%helium%logger)
      98          10 :             helium_env(k)%helium%logger => cp_get_default_logger()
      99             : 
     100             :             ! create iteration level
     101             :             ! Although the Helium MC is not really an md it is most of the times
     102             :             ! embedded in the pint code so the same step counter applies.
     103          10 :             IF (k .EQ. 1) THEN
     104          10 :                CALL cp_add_iter_level(helium_env(k)%helium%logger%iter_info, "PINT")
     105             :                CALL cp_iterate(helium_env(k)%helium%logger%iter_info, &
     106          10 :                                iter_nr=helium_env(k)%helium%first_step)
     107             :             END IF
     108          10 :             tot_steps = helium_env(k)%helium%first_step
     109          10 :             num_steps = helium_env(k)%helium%num_steps
     110             : 
     111             :             ! set the properties accumulated over the whole MC process to 0
     112          40 :             helium_env(k)%helium%proarea%accu(:) = 0.0_dp
     113          40 :             helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
     114          40 :             helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
     115          40 :             helium_env(k)%helium%mominer%accu(:) = 0.0_dp
     116          10 :             IF (helium_env(k)%helium%rho_present) THEN
     117        4222 :                helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
     118             :             END IF
     119          20 :             IF (helium_env(k)%helium%rdf_present) THEN
     120        1002 :                helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
     121             :             END IF
     122             :          END DO
     123             :       END IF
     124             : 
     125             :       ! Distribute steps for processors without helium_env
     126          10 :       CALL logger%para_env%bcast(num_steps)
     127          10 :       CALL logger%para_env%bcast(tot_steps)
     128             : 
     129          48 :       DO step = 1, num_steps
     130             : 
     131          38 :          tot_steps = tot_steps + 1
     132             : 
     133          38 :          IF (ASSOCIATED(helium_env)) THEN
     134          76 :             DO k = 1, SIZE(helium_env)
     135          38 :                IF (k .EQ. 1) THEN
     136             :                   CALL cp_iterate(helium_env(k)%helium%logger%iter_info, &
     137             :                                   last=(step == helium_env(k)%helium%num_steps), &
     138          38 :                                   iter_nr=tot_steps)
     139             :                END IF
     140          76 :                helium_env(k)%helium%current_step = tot_steps
     141             :             END DO
     142             :          END IF
     143             : 
     144          38 :          CALL helium_step(helium_env, pint_env)
     145             : 
     146             :          ! call write_restart here to avoid interference with PINT write_restart
     147          38 :          IF (ASSOCIATED(helium_env)) THEN
     148          38 :             CALL write_restart(root_section=helium_env(1)%helium%input, helium_env=helium_env)
     149             :          END IF
     150             : 
     151             :          ! exit from the main loop if soft exit has been requested
     152          38 :          CALL external_control(should_stop, "PINT", globenv=globenv)
     153          48 :          IF (should_stop) EXIT
     154             : 
     155             :       END DO
     156             : 
     157             :       ! remove iteration level
     158          10 :       IF (ASSOCIATED(helium_env)) THEN
     159          10 :          CALL cp_rm_iter_level(helium_env(1)%helium%logger%iter_info, "PINT")
     160             :       END IF
     161             : 
     162          10 :       RETURN
     163         320 :    END SUBROUTINE helium_do_run
     164             : 
     165             : ! ***************************************************************************
     166             : !> \brief  Sample the helium phase space
     167             : !> \param helium_env ...
     168             : !> \param pint_env ...
     169             : !> \date   2009-10-27
     170             : !> \par    History
     171             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     172             : !> \author Lukasz Walewski
     173             : !> \note   Samples helium variable space according to multilevel Metropolis
     174             : !>         MC scheme, calculates the forces exerted by helium solvent on the
     175             : !>         solute and stores them in helium%force_avrg array. The forces are
     176             : !>         averaged over outer MC loop.
     177             : !> \note   The implicit assumption is that we simulate solute _with_ solvent
     178             : !>         most of the time, so for performance reasons I do not check that
     179             : !>         everywhere I should. This leads to some redundancy in the case of
     180             : !>         helium-only calculations.
     181             : ! **************************************************************************************************
     182         117 :    SUBROUTINE helium_sample(helium_env, pint_env)
     183             : 
     184             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     185             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     186             : 
     187             :       CHARACTER(len=default_string_length)               :: msg_str
     188             :       INTEGER                                            :: i, irot, iweight, k, nslices, nsteps, &
     189             :                                                             num_env, offset, sel_mp_source
     190             :       REAL(KIND=dp)                                      :: inv_num_env, inv_xn, rnd, rtmp, rweight
     191         117 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work_2d
     192             :       TYPE(cp_logger_type), POINTER                      :: logger
     193             : 
     194         117 :       NULLIFY (logger)
     195         117 :       IF (SIZE(helium_env) < 1) THEN
     196           0 :          logger => cp_get_default_logger()
     197             :       ELSE
     198         117 :          CALL cp_logger_create(logger, para_env=helium_env(1)%comm, template_logger=cp_get_default_logger())
     199         117 :          CALL cp_add_default_logger(logger)
     200             :       END IF
     201             : 
     202         270 :       DO k = 1, SIZE(helium_env)
     203             : 
     204         153 :          IF (helium_env(k)%helium%solute_present) THEN
     205        7066 :             helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
     206         109 :             helium_env(k)%helium%current_step = pint_env%iter
     207         436 :             helium_env(k)%helium%origin = pint_com_pos(pint_env)
     208             :          END IF
     209             : 
     210        1683 :          helium_env(k)%helium%energy_avrg(:) = 0.0_dp
     211        3369 :          helium_env(k)%helium%plength_avrg(:) = 0.0_dp
     212        2178 :          helium_env(k)%helium%num_accepted(:, :) = 0.0_dp
     213             : 
     214             :          ! helium parallelization: each processor gets different RN stream and
     215             :          ! runs independent helium simulation, the properties and forces are
     216             :          ! averaged over parallel helium environments once per step.
     217         153 :          inv_xn = 0.0_dp
     218         180 :          SELECT CASE (helium_env(k)%helium%sampling_method)
     219             : 
     220             :          CASE (helium_sampling_worm)
     221             : 
     222          27 :             CALL helium_sample_worm(helium_env(k)%helium, pint_env)
     223             : 
     224          27 :             inv_xn = 1.0_dp/REAL(helium_env(k)%helium%worm_nstat, dp)
     225             : 
     226             :          CASE (helium_sampling_ceperley)
     227             :             ! helium sampling (outer MC loop)
     228        7146 :             DO irot = 1, helium_env(k)%helium%iter_rot
     229             : 
     230             :                ! rotate helium beads in imaginary time at random, store current
     231             :                ! 'rotation state' in helium%relrot which is within (0, helium%beads-1)
     232             :                ! (this is needed to sample different fragments of the permutation
     233             :                ! paths in try_permutations)
     234        7020 :                rnd = helium_env(k)%helium%rng_stream_uniform%next()
     235        7020 :                nslices = INT(rnd*helium_env(k)%helium%beads)
     236        7020 :                CALL helium_rotate(helium_env(k)%helium, nslices)
     237             : 
     238        7020 :                CALL helium_try_permutations(helium_env(k)%helium, pint_env)
     239             : 
     240             :                ! calculate & accumulate instantaneous properties for averaging
     241        7020 :                IF (helium_env(k)%helium%solute_present) THEN
     242        3620 :                   IF (helium_env(k)%helium%get_helium_forces == helium_forces_average) THEN
     243        2620 :                      CALL helium_solute_e_f(pint_env, helium_env(k)%helium, rtmp)
     244             :                      helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_avrg(:, :) + &
     245      120520 :                                                              helium_env(k)%helium%force_inst(:, :)
     246             :                   END IF
     247             :                END IF
     248        7020 :                CALL helium_calc_energy(helium_env(k)%helium, pint_env)
     249             : 
     250       77220 :                helium_env(k)%helium%energy_avrg(:) = helium_env(k)%helium%energy_avrg(:) + helium_env(k)%helium%energy_inst(:)
     251        7020 :                CALL helium_calc_plength(helium_env(k)%helium)
     252      220446 :                helium_env(k)%helium%plength_avrg(:) = helium_env(k)%helium%plength_avrg(:) + helium_env(k)%helium%plength_inst(:)
     253             : 
     254             :                ! instantaneous force output according to HELIUM%PRINT%FORCES_INST
     255             :                ! Warning: file I/O here may cost A LOT of cpu time!
     256             :                ! switched off here to save cpu
     257             :                !CALL helium_print_force_inst( helium_env(k)%helium, error )
     258             : 
     259             :             END DO ! outer MC loop
     260             : 
     261             :             ! If only the last forces are taken, calculate them for the last outer MC loop configuration
     262         126 :             IF ((helium_env(k)%helium%solute_present) .AND. (helium_env(k)%helium%get_helium_forces == helium_forces_last)) THEN
     263          10 :                CALL helium_solute_e_f(pint_env, helium_env(k)%helium, rtmp)
     264         460 :                helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_inst(:, :)
     265             :             END IF
     266             : 
     267             :             ! restore the original alignment of beads in imaginary time
     268             :             ! (this is useful to make a single bead's position easy to follow
     269             :             ! in the trajectory, otherwise its index would change every step)
     270         126 :             CALL helium_rotate(helium_env(k)%helium, -helium_env(k)%helium%relrot)
     271             : 
     272         126 :             inv_xn = 1.0_dp/REAL(helium_env(k)%helium%iter_rot, dp)
     273             : 
     274             :          CASE DEFAULT
     275           0 :             WRITE (msg_str, *) helium_env(k)%helium%sampling_method
     276           0 :             msg_str = "Sampling method ("//TRIM(ADJUSTL(msg_str))//") not supported."
     277         153 :             CPABORT(msg_str)
     278             : 
     279             :          END SELECT
     280             : 
     281             :          ! actually average the properties over the outer MC loop
     282         153 :          IF ((helium_env(k)%helium%solute_present) .AND. (helium_env(k)%helium%get_helium_forces == helium_forces_average)) THEN
     283        3772 :             helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_avrg(:, :)*inv_xn
     284             :          END IF
     285        1683 :          helium_env(k)%helium%energy_avrg(:) = helium_env(k)%helium%energy_avrg(:)*inv_xn
     286        3369 :          helium_env(k)%helium%plength_avrg(:) = helium_env(k)%helium%plength_avrg(:)*inv_xn
     287             : 
     288             :          ! instantaneous properties
     289         612 :          helium_env(k)%helium%proarea%inst(:) = helium_total_projected_area(helium_env(k)%helium)
     290         612 :          helium_env(k)%helium%prarea2%inst(:) = helium_env(k)%helium%proarea%inst(:)*helium_env(k)%helium%proarea%inst(:)
     291         612 :          helium_env(k)%helium%wnumber%inst(:) = helium_total_winding_number(helium_env(k)%helium)
     292         612 :          helium_env(k)%helium%wnmber2%inst(:) = helium_env(k)%helium%wnumber%inst(:)*helium_env(k)%helium%wnumber%inst(:)
     293         612 :          helium_env(k)%helium%mominer%inst(:) = helium_total_moment_of_inertia(helium_env(k)%helium)
     294             : 
     295             :          ! properties accumulated over the whole MC process
     296         612 :          helium_env(k)%helium%proarea%accu(:) = helium_env(k)%helium%proarea%accu(:) + helium_env(k)%helium%proarea%inst(:)
     297         612 :          helium_env(k)%helium%prarea2%accu(:) = helium_env(k)%helium%prarea2%accu(:) + helium_env(k)%helium%prarea2%inst(:)
     298         612 :          helium_env(k)%helium%wnmber2%accu(:) = helium_env(k)%helium%wnmber2%accu(:) + helium_env(k)%helium%wnmber2%inst(:)
     299         612 :          helium_env(k)%helium%mominer%accu(:) = helium_env(k)%helium%mominer%accu(:) + helium_env(k)%helium%mominer%inst(:)
     300         153 :          IF (helium_env(k)%helium%rho_present) THEN
     301          12 :             CALL helium_calc_rho(helium_env(k)%helium)
     302             :             helium_env(k)%helium%rho_accu(:, :, :, :) = helium_env(k)%helium%rho_accu(:, :, :, :) + &
     303       25332 :                                                         helium_env(k)%helium%rho_inst(:, :, :, :)
     304             :          END IF
     305         153 :          IF (helium_env(k)%helium%rdf_present) THEN
     306          12 :             CALL helium_set_rdf_coord_system(helium_env(k)%helium, pint_env)
     307          12 :             CALL helium_calc_rdf(helium_env(k)%helium)
     308        6012 :             helium_env(k)%helium%rdf_accu(:, :) = helium_env(k)%helium%rdf_accu(:, :) + helium_env(k)%helium%rdf_inst(:, :)
     309             :          END IF
     310             : 
     311             :          ! running averages (restart-aware)
     312         153 :          nsteps = helium_env(k)%helium%current_step - helium_env(k)%helium%first_step
     313         153 :          iweight = helium_env(k)%helium%averages_iweight
     314         153 :          rweight = REAL(iweight, dp)
     315         153 :          rtmp = 1.0_dp/(REAL(MAX(1, nsteps + iweight), dp))
     316             :          helium_env(k)%helium%proarea%ravr(:) = (helium_env(k)%helium%proarea%accu(:) + &
     317         612 :                                                  rweight*helium_env(k)%helium%proarea%rstr(:))*rtmp
     318             :          helium_env(k)%helium%prarea2%ravr(:) = (helium_env(k)%helium%prarea2%accu(:) + &
     319         612 :                                                  rweight*helium_env(k)%helium%prarea2%rstr(:))*rtmp
     320             :          helium_env(k)%helium%wnmber2%ravr(:) = (helium_env(k)%helium%wnmber2%accu(:) + &
     321         612 :                                                  rweight*helium_env(k)%helium%wnmber2%rstr(:))*rtmp
     322             :          helium_env(k)%helium%mominer%ravr(:) = (helium_env(k)%helium%mominer%accu(:) + &
     323         729 :                                                  rweight*helium_env(k)%helium%mominer%rstr(:))*rtmp
     324             : 
     325             :       END DO
     326             : 
     327             :       ! average over helium environments sitting at different processors
     328             : !TODO the following involves message passing, maybe move it to i/o routines?
     329         117 :       num_env = helium_env(1)%helium%num_env
     330         117 :       inv_num_env = 1.0_dp/REAL(num_env, dp)
     331             : 
     332             :       !energy_avrg:
     333         153 :       DO k = 2, SIZE(helium_env)
     334             :          helium_env(1)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:) + &
     335         513 :                                                helium_env(k)%helium%energy_avrg(:)
     336             :       END DO
     337         117 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%energy_avrg(:))
     338        1287 :       helium_env(1)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:)*inv_num_env
     339         153 :       DO k = 2, SIZE(helium_env)
     340         513 :          helium_env(k)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:)
     341             :       END DO
     342             : 
     343             :       !plength_avrg:
     344         153 :       DO k = 2, SIZE(helium_env)
     345             :          helium_env(1)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:) + &
     346         333 :                                                 helium_env(k)%helium%plength_avrg(:)
     347             :       END DO
     348        6189 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%plength_avrg(:))
     349        3153 :       helium_env(1)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:)*inv_num_env
     350         153 :       DO k = 2, SIZE(helium_env)
     351         333 :          helium_env(k)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:)
     352             :       END DO
     353             : 
     354             :       !num_accepted:
     355         153 :       DO k = 2, SIZE(helium_env)
     356             :          helium_env(1)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :) + &
     357         369 :                                                    helium_env(k)%helium%num_accepted(:, :)
     358             :       END DO
     359        3735 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%num_accepted(:, :))
     360        1926 :       helium_env(1)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :)*inv_num_env
     361         153 :       DO k = 2, SIZE(helium_env)
     362         369 :          helium_env(k)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :)
     363             :       END DO
     364             : 
     365             :       !force_avrg:
     366         117 :       IF (helium_env(1)%helium%solute_present) THEN
     367          73 :          IF (helium_env(1)%helium%get_helium_forces == helium_forces_average) THEN
     368          82 :             DO k = 2, SIZE(helium_env)
     369             :                helium_env(1)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :) + &
     370        1702 :                                                        helium_env(k)%helium%force_avrg(:, :)
     371             :             END DO
     372        4186 :             CALL helium_env(1)%comm%sum(helium_env(1)%helium%force_avrg(:, :))
     373        2116 :             helium_env(1)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :)*inv_num_env
     374          82 :             DO k = 2, SIZE(helium_env)
     375        1702 :                helium_env(k)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :)
     376             :             END DO
     377          27 :          ELSE IF (helium_env(1)%helium%get_helium_forces == helium_forces_last) THEN
     378          27 :             IF (logger%para_env%is_source()) THEN
     379             :                sel_mp_source = INT(helium_env(1)%helium%rng_stream_uniform%next() &
     380          16 :                                    *REAL(helium_env(1)%helium%num_env, dp))
     381             :             END IF
     382          27 :             CALL helium_env(1)%comm%bcast(sel_mp_source, logger%para_env%source)
     383             : 
     384          27 :             offset = 0
     385          38 :             DO i = 1, logger%para_env%mepos
     386          38 :                offset = offset + helium_env(1)%env_all(i)
     387             :             END DO
     388             :             ALLOCATE (work_2d(SIZE(helium_env(1)%helium%force_avrg, 1), &
     389         108 :                               SIZE(helium_env(1)%helium%force_avrg, 2)))
     390        3294 :             work_2d(:, :) = 0.0_dp
     391          27 :             IF (sel_mp_source .GE. offset .AND. sel_mp_source .LT. offset + SIZE(helium_env)) THEN
     392        2032 :                work_2d(:, :) = helium_env(sel_mp_source - offset + 1)%helium%force_avrg(:, :)
     393             :             END IF
     394          27 :             CALL helium_env(1)%comm%sum(work_2d(:, :))
     395          54 :             DO k = 1, SIZE(helium_env)
     396        3321 :                helium_env(k)%helium%force_avrg(:, :) = work_2d(:, :)
     397             :             END DO
     398          27 :             DEALLOCATE (work_2d)
     399             :          END IF
     400             :       END IF
     401             : 
     402         117 :       IF (SIZE(helium_env) > 0) THEN
     403         117 :          CALL cp_rm_default_logger()
     404         117 :          CALL cp_logger_release(logger)
     405             :       END IF
     406             : 
     407         117 :       RETURN
     408             :    END SUBROUTINE helium_sample
     409             : 
     410             : ! ***************************************************************************
     411             : !> \brief  Perform MC step for helium
     412             : !> \param helium_env ...
     413             : !> \param pint_env ...
     414             : !> \date   2009-06-12
     415             : !> \par    History
     416             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     417             : !> \author Lukasz Walewski
     418             : ! **************************************************************************************************
     419         110 :    SUBROUTINE helium_step(helium_env, pint_env)
     420             : 
     421             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     422             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     423             : 
     424             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_step'
     425             : 
     426             :       CHARACTER(len=default_string_length)               :: msgstr, stmp, time_unit
     427             :       INTEGER                                            :: handle, k
     428             :       REAL(KIND=dp)                                      :: time_start, time_stop, time_used
     429         110 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: DATA
     430             : 
     431         110 :       CALL timeset(routineN, handle)
     432         110 :       time_start = m_walltime()
     433             : 
     434         110 :       IF (ASSOCIATED(helium_env)) THEN
     435             :          ! perform the actual phase space sampling
     436         105 :          CALL helium_sample(helium_env, pint_env)
     437             : 
     438             :          ! print properties
     439         105 :          CALL helium_print_energy(helium_env)
     440         105 :          CALL helium_print_plength(helium_env)
     441         105 :          CALL helium_print_accepts(helium_env)
     442         105 :          CALL helium_print_force(helium_env)
     443         105 :          CALL helium_print_perm(helium_env)
     444         105 :          CALL helium_print_coordinates(helium_env)
     445         105 :          CALL helium_print_action(pint_env, helium_env)
     446             : 
     447         105 :          IF (helium_env(1)%helium%rdf_present) CALL helium_print_rdf(helium_env)
     448         105 :          IF (helium_env(1)%helium%rho_present) CALL helium_print_rho(helium_env)
     449             : 
     450         105 :          NULLIFY (DATA)
     451         315 :          ALLOCATE (DATA(3*SIZE(helium_env)))
     452             : 
     453         105 :          WRITE (stmp, *) helium_env(1)%helium%apref
     454         510 :          DATA(:) = 0.0_dp
     455         240 :          DO k = 1, SIZE(helium_env)
     456         645 :             DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%proarea%inst(:)
     457             :          END DO
     458             :          CALL helium_print_vector(helium_env, &
     459             :                                   "MOTION%PINT%HELIUM%PRINT%PROJECTED_AREA", &
     460             :                                   DATA, &
     461             :                                   angstrom*angstrom, &
     462             :                                   (/"A_x", "A_y", "A_z"/), &
     463             :                                   "prefactor = "//TRIM(ADJUSTL(stmp))//" [Angstrom^-2]", &
     464         420 :                                   "helium-pa")
     465             : 
     466         510 :          DATA(:) = 0.0_dp
     467         240 :          DO k = 1, SIZE(helium_env)
     468         645 :             DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%prarea2%ravr(:)
     469             :          END DO
     470             :          CALL helium_print_vector(helium_env, &
     471             :                                   "MOTION%PINT%HELIUM%PRINT%PROJECTED_AREA_2_AVG", &
     472             :                                   DATA, &
     473             :                                   angstrom*angstrom*angstrom*angstrom, &
     474             :                                   (/"<A_x^2>", "<A_y^2>", "<A_z^2>"/), &
     475             :                                   "prefactor = "//TRIM(ADJUSTL(stmp))//" [Angstrom^-2]", &
     476             :                                   "helium-pa-avg", &
     477             :                                   "REWIND", &
     478         420 :                                   .TRUE.)
     479             : 
     480         510 :          DATA(:) = 0.0_dp
     481         240 :          DO k = 1, SIZE(helium_env)
     482         645 :             DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%mominer%inst(:)
     483             :          END DO
     484             :          CALL helium_print_vector(helium_env, &
     485             :                                   "MOTION%PINT%HELIUM%PRINT%MOMENT_OF_INERTIA", &
     486             :                                   DATA, &
     487             :                                   angstrom*angstrom, &
     488             :                                   (/"I_x/m", "I_y/m", "I_z/m"/), &
     489             :                                   "prefactor = "//TRIM(ADJUSTL(stmp))//" [Angstrom^-2]", &
     490         420 :                                   "helium-mi")
     491             : 
     492         510 :          DATA(:) = 0.0_dp
     493         240 :          DO k = 1, SIZE(helium_env)
     494         645 :             DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%mominer%ravr
     495             :          END DO
     496             :          CALL helium_print_vector(helium_env, &
     497             :                                   "MOTION%PINT%HELIUM%PRINT%MOMENT_OF_INERTIA_AVG", &
     498             :                                   DATA, &
     499             :                                   angstrom*angstrom, &
     500             :                                   (/"<I_x/m>", "<I_y/m>", "<I_z/m>"/), &
     501             :                                   "prefactor = "//TRIM(ADJUSTL(stmp))//" [Angstrom^-2]", &
     502             :                                   "helium-mi-avg", &
     503             :                                   "REWIND", &
     504         420 :                                   .TRUE.)
     505             : 
     506         510 :          DATA(:) = 0.0_dp
     507         240 :          DO k = 1, SIZE(helium_env)
     508         645 :             DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%wnumber%inst
     509             :          END DO
     510         105 :          WRITE (stmp, *) helium_env(1)%helium%wpref
     511             :          CALL helium_print_vector(helium_env, &
     512             :                                   "MOTION%PINT%HELIUM%PRINT%WINDING_NUMBER", &
     513             :                                   DATA, &
     514             :                                   angstrom, &
     515             :                                   (/"W_x", "W_y", "W_z"/), &
     516             :                                   "prefactor = "//TRIM(ADJUSTL(stmp))//" [Angstrom^-2]", &
     517         420 :                                   "helium-wn")
     518             : 
     519         510 :          DATA(:) = 0.0_dp
     520         240 :          DO k = 1, SIZE(helium_env)
     521         645 :             DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%wnmber2%ravr
     522             :          END DO
     523             :          CALL helium_print_vector(helium_env, &
     524             :                                   "MOTION%PINT%HELIUM%PRINT%WINDING_NUMBER_2_AVG", &
     525             :                                   DATA, &
     526             :                                   angstrom*angstrom, &
     527             :                                   (/"<W_x^2>", "<W_y^2>", "<W_z^2>"/), &
     528             :                                   "prefactor = "//TRIM(ADJUSTL(stmp))//" [Angstrom^-2]", &
     529             :                                   "helium-wn-avg", &
     530             :                                   "REWIND", &
     531         420 :                                   .TRUE.)
     532         105 :          DEALLOCATE (DATA)
     533             : 
     534         105 :          time_stop = m_walltime()
     535         105 :          time_used = time_stop - time_start
     536         105 :          time_unit = "sec"
     537         105 :          IF (time_used .GE. 60.0_dp) THEN
     538           0 :             time_used = time_used/60.0_dp
     539           0 :             time_unit = "min"
     540           0 :             IF (time_used .GE. 60.0_dp) THEN
     541           0 :                time_used = time_used/60.0_dp
     542           0 :                time_unit = "hours"
     543             :             END IF
     544             :          END IF
     545         105 :          msgstr = "MC step"
     546         105 :          stmp = ""
     547         105 :          WRITE (stmp, *) helium_env(1)%helium%current_step
     548         105 :          msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))//" of"
     549         105 :          stmp = ""
     550         105 :          WRITE (stmp, *) helium_env(1)%helium%last_step
     551         105 :          msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))//" in"
     552         105 :          stmp = ""
     553         105 :          WRITE (stmp, '(F20.1)') time_used
     554         105 :          msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))
     555         105 :          msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(time_unit))//"."
     556         105 :          CALL helium_write_line(TRIM(msgstr))
     557             : 
     558             :          ! print out the total energy - for regtest evaluation
     559         105 :          stmp = ""
     560         105 :          WRITE (stmp, *) helium_env(1)%helium%energy_avrg(e_id_total)
     561         105 :          msgstr = "Total energy = "//TRIM(ADJUSTL(stmp))
     562         105 :          CALL helium_write_line(TRIM(msgstr))
     563             :       END IF
     564             : 
     565         110 :       CALL timestop(handle)
     566             : 
     567         110 :       RETURN
     568         110 :    END SUBROUTINE helium_step
     569             : 
     570             : ! ***************************************************************************
     571             : !> \brief  ...
     572             : !> \param helium ...
     573             : !> \param  pint_env    path integral environment
     574             : !> \par    History
     575             : !>        2010-06-17 added different distributions for m-sampling [lwalewski]
     576             : !>        2010-06-17 ratio for m-value added (m-sampling related) [lwalewski]
     577             : !> \author hforbert
     578             : ! **************************************************************************************************
     579        7020 :    SUBROUTINE helium_try_permutations(helium, pint_env)
     580             :       TYPE(helium_solvent_type), POINTER                 :: helium
     581             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     582             : 
     583             :       CHARACTER(len=default_string_length)               :: err_str, stmp
     584             :       INTEGER                                            :: cyclen, ni
     585             :       LOGICAL                                            :: accepted, selected
     586             :       REAL(KIND=dp)                                      :: r, rnd, x, y, z
     587             : 
     588        7020 :       IF (helium%maxcycle > 1) CALL helium_update_transition_matrix(helium)
     589             : 
     590             :       ! consecutive calls to helium_slice_metro_cyclic require that
     591             :       ! ALL(helium%pos.EQ.helium%work) - see a comment there,
     592             :       ! besides there is a magic regarding helium%permutation
     593             :       ! you have been warned
     594    17717940 :       helium%work(:, :, :) = helium%pos(:, :, :)
     595             : 
     596             :       ! the inner MC loop (without rotation in imaginary time)
     597     1348020 :       DO ni = 1, helium%iter_norot
     598             : 
     599             :          ! set the probability threshold for m_value: 1/(1+(m-1)/helium%m_ratio)
     600     1341000 :          r = 1.0_dp/(1.0_dp + (helium%maxcycle - 1)/helium%m_ratio)
     601             : 
     602             :          ! draw permutation length for this trial from the distribution of choice
     603             :          !
     604     1341000 :          SELECT CASE (helium%m_dist_type)
     605             : 
     606             :          CASE (helium_mdist_singlev)
     607           0 :             x = helium%rng_stream_uniform%next()
     608           0 :             IF (x .LT. r) THEN
     609           0 :                cyclen = 1
     610             :             ELSE
     611           0 :                cyclen = helium%m_value
     612             :             END IF
     613             : 
     614             :          CASE (helium_mdist_uniform)
     615     1341000 :             x = helium%rng_stream_uniform%next()
     616     1341000 :             IF (x .LT. r) THEN
     617      350513 :                cyclen = helium%m_value
     618             :             ELSE
     619             :                DO
     620     1320591 :                   x = helium%rng_stream_uniform%next()
     621     1320591 :                   cyclen = INT(helium%maxcycle*x) + 1
     622     1320591 :                   IF (cyclen .NE. helium%m_value) EXIT
     623             :                END DO
     624             :             END IF
     625             : 
     626             :          CASE (helium_mdist_linear)
     627           0 :             x = helium%rng_stream_uniform%next()
     628           0 :             IF (x .LT. r) THEN
     629           0 :                cyclen = helium%m_value
     630             :             ELSE
     631             :                DO
     632           0 :                   x = helium%rng_stream_uniform%next()
     633           0 :                   y = SQRT(2.0_dp*x)
     634           0 :                   cyclen = INT(helium%maxcycle*y/SQRT(2.0_dp)) + 1
     635           0 :                   IF (cyclen .NE. helium%m_value) EXIT
     636             :                END DO
     637             :             END IF
     638             : 
     639             :          CASE (helium_mdist_quadratic)
     640           0 :             x = helium%rng_stream_uniform%next()
     641           0 :             IF (x .LT. r) THEN
     642           0 :                cyclen = helium%m_value
     643             :             ELSE
     644             :                DO
     645           0 :                   x = helium%rng_stream_uniform%next()
     646           0 :                   y = (3.0_dp*x)**(1.0_dp/3.0_dp)
     647           0 :                   z = 3.0_dp**(1.0_dp/3.0_dp)
     648           0 :                   cyclen = INT(helium%maxcycle*y/z) + 1
     649           0 :                   IF (cyclen .NE. helium%m_value) EXIT
     650             :                END DO
     651             :             END IF
     652             : 
     653             :          CASE (helium_mdist_exponential)
     654           0 :             x = helium%rng_stream_uniform%next()
     655           0 :             IF (x .LT. r) THEN
     656           0 :                cyclen = helium%m_value
     657             :             ELSE
     658             :                DO
     659             :                   DO
     660           0 :                      x = helium%rng_stream_uniform%next()
     661           0 :                      IF (x .GE. 0.01_dp) EXIT
     662             :                   END DO
     663           0 :                   z = -LOG(0.01_dp)
     664           0 :                   y = LOG(x)/z + 1.0_dp; 
     665           0 :                   cyclen = INT(helium%maxcycle*y) + 1
     666           0 :                   IF (cyclen .NE. helium%m_value) EXIT
     667             :                END DO
     668             :             END IF
     669             : 
     670             :          CASE (helium_mdist_gaussian)
     671           0 :             x = helium%rng_stream_uniform%next()
     672           0 :             IF (x .LT. r) THEN
     673           0 :                cyclen = 1
     674             :             ELSE
     675             :                DO
     676           0 :                   x = helium%rng_stream_gaussian%next()
     677           0 :                   cyclen = INT(x*0.75_dp + helium%m_value - 0.5_dp) + 1
     678           0 :                   IF (cyclen .NE. 1) EXIT
     679             :                END DO
     680             :             END IF
     681             : 
     682             :          CASE DEFAULT
     683           0 :             WRITE (stmp, *) helium%m_dist_type
     684           0 :             err_str = "M distribution type unknown ("//TRIM(ADJUSTL(stmp))//")"
     685           0 :             CPABORT(err_str)
     686     1341000 :             cyclen = 1
     687             : 
     688             :          END SELECT
     689             : 
     690     1341000 :          IF (cyclen < 1) cyclen = 1
     691     1341000 :          IF (cyclen > helium%maxcycle) cyclen = helium%maxcycle
     692     1341000 :          helium%num_accepted(1, cyclen) = helium%num_accepted(1, cyclen) + 1
     693             : 
     694             :          ! check, if permutation of this length can be constructed
     695     1341000 :          IF (cyclen == 1) THEN
     696      350513 :             rnd = helium%rng_stream_uniform%next()
     697      350513 :             helium%ptable(1) = 1 + INT(rnd*helium%atoms)
     698      350513 :             helium%ptable(2) = -1
     699      350513 :             helium%pweight = 0.0_dp
     700             :             selected = .TRUE.
     701             :          ELSE
     702      990487 :             selected = helium_select_permutation(helium, cyclen)
     703             :          END IF
     704             : 
     705      990487 :          IF (selected) THEN
     706             :             ! permutation was successfully selected - actually sample this permutation
     707      350548 :             accepted = helium_slice_metro_cyclic(helium, pint_env, cyclen)
     708             :          ELSE
     709             :             accepted = .FALSE.
     710             :          END IF
     711             : 
     712             : !if (any(helium%pos .ne. helium%work)) then
     713             : !  print *, ""
     714             : !  print *, "pos and work are different!"
     715             : !  print *, ""
     716             : !end if
     717             : 
     718      357568 :          IF (accepted) THEN
     719             :             ! positions changed
     720      128676 :             IF (helium%solute_present) THEN
     721             :                ! use COM of the solute, but it did not move here so do nothing to save cpu
     722             :             ELSE
     723       60356 :                IF (helium%periodic) THEN
     724             :                   ! use center of the cell, but it did not move here so do nothing to save cpu
     725             :                ELSE
     726             :                   ! update center of mass
     727           0 :                   helium%center(:) = helium_com(helium)
     728             :                END IF
     729             :             END IF
     730             :          END IF
     731             : 
     732             :       END DO
     733             : 
     734        7020 :       RETURN
     735             :    END SUBROUTINE helium_try_permutations
     736             : 
     737             : ! ***************************************************************************
     738             : !> \brief  Sample path variables for permutation of length <cyclen>
     739             : !> \param helium ...
     740             : !> \param pint_env ...
     741             : !> \param cyclen ...
     742             : !> \return ...
     743             : !> \author hforbert
     744             : ! **************************************************************************************************
     745      350548 :    FUNCTION helium_slice_metro_cyclic(helium, pint_env, cyclen) RESULT(res)
     746             :       TYPE(helium_solvent_type), POINTER                 :: helium
     747             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     748             :       INTEGER, INTENT(IN)                                :: cyclen
     749             :       LOGICAL                                            :: res
     750             : 
     751             :       CHARACTER(len=default_string_length)               :: err_str, stmp
     752             :       INTEGER                                            :: c, ia, ib, ic, ifix, j, k, l, level, &
     753             :                                                             pk1, pk2, stride
     754      350548 :       INTEGER, DIMENSION(:), POINTER                     :: p, perm
     755             :       LOGICAL                                            :: nperiodic, should_reject
     756             :       REAL(KIND=dp)                                      :: cell_size, ds, dtk, e1, e2, pds, &
     757             :                                                             prev_ds, r, rtmp, rtmpo, sigma, x
     758      350548 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     759             :       REAL(KIND=dp), DIMENSION(3)                        :: bis, biso, new_com, rm1, rm2, tmp1, tmp2
     760      350548 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos, work
     761             :       TYPE(spline_data_type), POINTER                    :: u0
     762             : 
     763             : ! trial permutation implicit in p
     764             : ! since we (momentarily) only use cyclic permutations:
     765             : ! cyclen = 1 : no permutation, sample p[0] anew
     766             : ! cyclen = 2 : p[0] -> p[1] -> p[0]
     767             : ! cyclen = 3 : p[0] -> p[1] -> p[2] -> p[0]
     768             : ! cyclen = 4 : p[0] -> p[1] -> p[2] -> p[3] -> p[0]
     769             : 
     770      350548 :       p => helium%ptable
     771      350548 :       prev_ds = helium%pweight
     772             : 
     773      350548 :       helium%num_accepted(2, cyclen) = helium%num_accepted(2, cyclen) + 1
     774      350548 :       level = 1
     775      350548 :       res = .FALSE.
     776             : 
     777             :       ! nothing to be done
     778      350548 :       IF (helium%bisection == 0) RETURN
     779             : 
     780             :       ! On entry helium_slice_metro_cyclic ASSUMES that work is a copy of pos
     781             :       ! and constructs the trial move there. If the move is accepted, the
     782             :       ! changed parts are copied to pos, but if it fails, only the CHANGED parts
     783             :       ! of work are overwritten by the corresponding parts of pos. So on exit work
     784             :       ! will AGAIN be a copy of pos (either way accept/reject). This is done here
     785             :       ! so we do not have to copy around the whole pos array in the caller, and
     786             :       ! the caller does not need to know which parts of work might have changed.
     787      350548 :       pos => helium%pos
     788      350548 :       work => helium%work
     789      350548 :       perm => helium%permutation
     790      350548 :       u0 => helium%u0
     791      350548 :       cell_size = (0.5_dp*helium%cell_size)**2
     792      350548 :       nperiodic = .NOT. helium%periodic
     793             : 
     794      350548 :       pds = prev_ds
     795      350548 :       ifix = helium%beads - helium%bisection + 1
     796             : 
     797             :       ! sanity checks
     798             :       !
     799      350548 :       IF (ifix < 1) THEN
     800           0 :          WRITE (stmp, *) ifix
     801           0 :          err_str = "ifix<1 test failed (ifix="//TRIM(ADJUSTL(stmp))//")"
     802           0 :          CPABORT(err_str)
     803             :       END IF
     804             :       !
     805      350548 :       j = 1
     806      350548 :       k = helium%bisection
     807     1051644 :       DO
     808     1402192 :          IF (k < 2) EXIT
     809     1051644 :          j = j*2
     810     1051644 :          k = k/2
     811             :       END DO
     812             :       !
     813      350548 :       IF (j /= helium%bisection) THEN
     814           0 :          WRITE (stmp, *) helium%bisection
     815           0 :          err_str = "helium%bisection not a power of 2 (helium%bisection="//TRIM(ADJUSTL(stmp))//")"
     816           0 :          CPABORT(err_str)
     817             :       END IF
     818             :       !
     819      350548 :       IF (helium%bisection < 2) THEN
     820           0 :          WRITE (stmp, *) helium%bisection
     821           0 :          err_str = "helium%bisection less than 2 (helium%bisection="//TRIM(ADJUSTL(stmp))//")"
     822           0 :          CPABORT(err_str)
     823             :       END IF
     824             : 
     825      350548 :       stride = helium%bisection
     826      399882 :       DO
     827      750430 :          IF (stride <= 2) EXIT
     828             :          ! calc new trial positions
     829             :          ! trial action: 0.5*stride*endpointapprox
     830      585672 :          sigma = SQRT(0.5_dp*helium%hb2m*(stride/2)*helium%tau)
     831      585672 :          dtk = 0.0_dp
     832      585672 :          ds = 0.0_dp
     833             : 
     834      585672 :          j = ifix + stride/2
     835      235124 :          DO
     836      820796 :             IF (j > helium%beads - stride/2) EXIT
     837      235124 :             pk1 = j - stride/2
     838      235124 :             pk2 = j + stride/2
     839             :             ! calculate log(T(s)):
     840      470250 :             DO k = 1, cyclen
     841      235126 :                CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, p(k), pk2), bis)
     842      940504 :                tmp1(:) = bis(:) - pos(:, p(k), j)
     843      235126 :                CALL helium_pbc(helium, tmp1)
     844      940504 :                tmp1(:) = tmp1(:)/sigma
     845      470250 :                dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     846             :             END DO
     847             :             ! calculate log(T(sprime)) and sprime itself
     848      470250 :             DO k = 1, cyclen
     849      235126 :                CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
     850      940504 :                DO c = 1, 3
     851      705378 :                   x = helium%rng_stream_gaussian%next(variance=1.0_dp)
     852      705378 :                   x = sigma*x
     853      705378 :                   tmp1(c) = tmp1(c) + x
     854      940504 :                   tmp2(c) = x
     855             :                END DO
     856      235126 :                CALL helium_pbc(helium, tmp1)
     857      235126 :                CALL helium_pbc(helium, tmp2)
     858      940504 :                work(:, p(k), j) = tmp1(:)
     859      940504 :                tmp2(:) = tmp2(:)/sigma
     860      470250 :                dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
     861             :             END DO
     862      235124 :             j = j + stride
     863             :          END DO
     864             : 
     865      585672 :          j = helium%beads - stride/2 + 1
     866      585672 :          pk1 = j - stride/2
     867     1171381 :          DO k = 1, cyclen
     868      585709 :             CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, perm(p(k)), 1), bis)
     869     2342836 :             tmp1(:) = bis(:) - pos(:, p(k), j)
     870      585709 :             CALL helium_pbc(helium, tmp1)
     871     2342836 :             tmp1(:) = tmp1(:)/sigma
     872     1171381 :             dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     873             :          END DO
     874     1171381 :          DO k = 1, cyclen
     875      585709 :             CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + MOD(k, cyclen))), 1), tmp1)
     876     2342836 :             DO c = 1, 3
     877     1757127 :                x = helium%rng_stream_gaussian%next(variance=1.0_dp)
     878     1757127 :                x = sigma*x
     879     1757127 :                tmp1(c) = tmp1(c) + x
     880     2342836 :                tmp2(c) = x
     881             :             END DO
     882      585709 :             CALL helium_pbc(helium, tmp1)
     883      585709 :             CALL helium_pbc(helium, tmp2)
     884     2342836 :             work(:, p(k), j) = tmp1(:)
     885     2342836 :             tmp2(:) = tmp2(:)/sigma
     886     1171381 :             dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
     887             :          END DO
     888             :          ! ok we got the new positions
     889             :          ! calculate action_k(s)-action_k(sprime)
     890      585672 :          x = 1.0_dp/(helium%tau*helium%hb2m*stride)
     891      585672 :          j = ifix
     892     1055920 :          DO
     893     1641592 :             IF (j > helium%beads - stride/2) EXIT
     894     1055920 :             pk1 = j + stride/2
     895     2111881 :             DO k = 1, cyclen
     896     4223844 :                tmp1(:) = pos(:, p(k), j) - pos(:, p(k), pk1)
     897     1055961 :                CALL helium_pbc(helium, tmp1)
     898     1055961 :                ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     899     4223844 :                tmp1(:) = work(:, p(k), j) - work(:, p(k), pk1)
     900     1055961 :                CALL helium_pbc(helium, tmp1)
     901     1055961 :                ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     902             :                ! interaction change
     903     1055961 :                IF (helium%solute_present) THEN
     904      543746 :                   CALL helium_bead_solute_e_f(pint_env, helium, p(k), pk1, energy=e1)
     905      543746 :                   CALL helium_bead_solute_e_f(pint_env, helium, p(k), pk1, work(:, p(k), pk1), e2)
     906      543746 :                   ds = ds + (stride/2)*(e1 - e2)*helium%tau
     907             :                END IF
     908    32647563 :                DO l = 1, helium%atoms
     909    32647563 :                   IF (l /= p(k)) THEN
     910   122142564 :                      tmp1(:) = pos(:, p(k), pk1) - pos(:, l, pk1)
     911    30535641 :                      CALL helium_pbc(helium, tmp1)
     912    30535641 :                      r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
     913    30535641 :                      IF ((r < cell_size) .OR. nperiodic) THEN
     914    28604255 :                         r = SQRT(r)
     915    28604255 :                         ds = ds + REAL(stride/2, dp)*helium_spline(u0, r)
     916             :                      END IF
     917   122142564 :                      tmp1(:) = work(:, p(k), pk1) - work(:, l, pk1)
     918    30535641 :                      CALL helium_pbc(helium, tmp1)
     919    30535641 :                      r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
     920    30535641 :                      IF ((r < cell_size) .OR. nperiodic) THEN
     921    28603810 :                         r = SQRT(r)
     922    28603810 :                         ds = ds - REAL(stride/2, dp)*helium_spline(u0, r)
     923             :                      END IF
     924             :                   END IF
     925             :                END DO
     926             :                ! counted p[k], p[m] twice. subtract those again
     927     2111881 :                IF (k < cyclen) THEN
     928          82 :                   DO l = k + 1, cyclen
     929         164 :                      tmp1(:) = pos(:, p(k), pk1) - pos(:, p(l), pk1)
     930          41 :                      CALL helium_pbc(helium, tmp1)
     931          41 :                      r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
     932          41 :                      IF ((r < cell_size) .OR. nperiodic) THEN
     933          41 :                         r = SQRT(r)
     934          41 :                         ds = ds - REAL(stride/2, dp)*helium_spline(u0, r)
     935             :                      END IF
     936         164 :                      tmp1(:) = work(:, p(k), pk1) - work(:, p(l), pk1)
     937          41 :                      CALL helium_pbc(helium, tmp1)
     938          41 :                      r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
     939          82 :                      IF ((r < cell_size) .OR. nperiodic) THEN
     940          41 :                         r = SQRT(r)
     941          41 :                         ds = ds + REAL(stride/2, dp)*helium_spline(u0, r)
     942             :                      END IF
     943             :                   END DO
     944             :                END IF
     945             :             END DO
     946     1055920 :             j = j + stride/2
     947             :          END DO
     948             :          ! last link
     949      585672 :          pk1 = helium%beads - stride/2 + 1
     950     1171381 :          DO k = 1, cyclen
     951     2342836 :             tmp1(:) = pos(:, p(k), pk1) - pos(:, perm(p(k)), 1)
     952      585709 :             CALL helium_pbc(helium, tmp1)
     953      585709 :             ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     954     2342836 :             tmp1(:) = work(:, p(k), pk1) - work(:, perm(p(1 + MOD(k, cyclen))), 1)
     955      585709 :             CALL helium_pbc(helium, tmp1)
     956     1171381 :             ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     957             :          END DO
     958             :          ! ok now accept or reject:
     959      585672 :          rtmp = helium%rng_stream_uniform%next()
     960             : !      IF ((dtk+ds-pds < 0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
     961      585672 :          IF (dtk + ds - pds < 0.0_dp) THEN
     962      376906 :             IF (EXP(dtk + ds - pds) < rtmp) THEN
     963     1672110 :                DO c = ifix, helium%beads
     964     3158710 :                   DO k = 1, cyclen
     965    11892520 :                      work(:, p(k), c) = pos(:, p(k), c)
     966             :                   END DO
     967             :                END DO
     968             :                RETURN
     969             :             END IF
     970             :          END IF
     971             :          ! accepted. go on to the next level
     972      399882 :          helium%num_accepted(level + 2, cyclen) = helium%num_accepted(level + 2, cyclen) + 1
     973      399882 :          level = level + 1
     974      399882 :          pds = ds
     975      399882 :          stride = stride/2
     976             :       END DO
     977             :       ! we are on the lowest level now
     978             :       ! calc new trial positions
     979             :       ! trial action: endpointapprox for T, full action otherwise
     980      164758 :       sigma = SQRT(0.5_dp*helium%hb2m*helium%tau)
     981      164758 :       dtk = 0.0_dp
     982      164758 :       ds = 0.0_dp
     983      659032 :       DO j = ifix + 1, helium%beads - 1, 2
     984      494274 :          pk1 = j - 1
     985      494274 :          pk2 = j + 1
     986             :          ! calculate log(T(s)):
     987      988548 :          DO k = 1, cyclen
     988      494274 :             CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, p(k), pk2), bis)
     989     1977096 :             tmp1(:) = bis(:) - pos(:, p(k), j)
     990      494274 :             CALL helium_pbc(helium, tmp1)
     991     1977096 :             tmp1(:) = tmp1(:)/sigma
     992      988548 :             dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
     993             :          END DO
     994             :          ! calculate log(T(sprime)) and sprime itself
     995     1153306 :          DO k = 1, cyclen
     996      494274 :             CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
     997     1977096 :             DO c = 1, 3
     998     1482822 :                x = helium%rng_stream_gaussian%next(variance=1.0_dp)
     999     1482822 :                x = sigma*x
    1000     1482822 :                tmp1(c) = tmp1(c) + x
    1001     1977096 :                tmp2(c) = x
    1002             :             END DO
    1003      494274 :             CALL helium_pbc(helium, tmp1)
    1004      494274 :             CALL helium_pbc(helium, tmp2)
    1005     1977096 :             work(:, p(k), j) = tmp1(:)
    1006     1977096 :             tmp2(:) = tmp2(:)/sigma
    1007      988548 :             dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
    1008             :          END DO
    1009             :       END DO
    1010      164758 :       j = helium%beads
    1011      164758 :       pk1 = j - 1
    1012      329516 :       DO k = 1, cyclen
    1013      164758 :          CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, perm(p(k)), 1), bis)
    1014      659032 :          tmp1(:) = bis(:) - pos(:, p(k), j)
    1015      164758 :          CALL helium_pbc(helium, tmp1)
    1016      659032 :          tmp1(:) = tmp1(:)/sigma
    1017      329516 :          dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
    1018             :       END DO
    1019      329516 :       DO k = 1, cyclen
    1020      164758 :          CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + MOD(k, cyclen))), 1), tmp1)
    1021      659032 :          DO c = 1, 3
    1022      494274 :             x = helium%rng_stream_gaussian%next(variance=1.0_dp)
    1023      494274 :             x = sigma*x
    1024      494274 :             tmp1(c) = tmp1(c) + x
    1025      659032 :             tmp2(c) = x
    1026             :          END DO
    1027      164758 :          CALL helium_pbc(helium, tmp1)
    1028      164758 :          CALL helium_pbc(helium, tmp2)
    1029      659032 :          work(:, p(k), j) = tmp1(:)
    1030      659032 :          tmp2 = tmp2/sigma
    1031      329516 :          dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
    1032             :       END DO
    1033             :       ! ok we got the new positions.
    1034             :       ! calculate action_k(s)-action_k(sprime)
    1035             :       ! interaction change
    1036             : !TODO interaction ONLY here? or some simplified 12-6 in the upper part?
    1037      164758 :       IF (helium%solute_present) THEN
    1038      772362 :          DO j = ifix, helium%beads
    1039     1458906 :             DO k = 1, cyclen
    1040      686544 :                CALL helium_bead_solute_e_f(pint_env, helium, p(k), j, energy=e1)
    1041      686544 :                CALL helium_bead_solute_e_f(pint_env, helium, p(k), j, work(:, p(k), j), e2)
    1042     2059632 :                ds = ds + (e1 - e2)*helium%tau
    1043             :             END DO
    1044             :          END DO
    1045             :       END IF
    1046      494274 :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
    1047      164758 :       x = 1.0_dp/(helium%tau*helium%hb2m*stride)
    1048     1318064 :       DO j = ifix, helium%beads - 1
    1049     1153306 :          pk1 = j + 1
    1050     2471370 :          DO k = 1, cyclen
    1051     4613224 :             tmp1(:) = pos(:, p(k), j) - pos(:, p(k), pk1)
    1052     1153306 :             CALL helium_pbc(helium, tmp1)
    1053     1153306 :             ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
    1054     4613224 :             tmp1(:) = work(:, p(k), j) - work(:, p(k), pk1)
    1055     1153306 :             CALL helium_pbc(helium, tmp1)
    1056     1153306 :             ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
    1057    34347327 :             DO l = 1, helium%atoms
    1058    34347327 :                IF (l /= p(k)) THEN
    1059   128162860 :                   rm1(:) = pos(:, p(k), j) - pos(:, l, j)
    1060   128162860 :                   rm2(:) = pos(:, p(k), pk1) - pos(:, l, pk1)
    1061    32040715 :                   ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
    1062   128162860 :                   rm1(:) = work(:, p(k), j) - work(:, l, j)
    1063   128162860 :                   rm2(:) = work(:, p(k), pk1) - work(:, l, pk1)
    1064    32040715 :                   ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
    1065             :                END IF
    1066             :             END DO
    1067             :             ! counted p[k], p[m] twice. subtract those again
    1068     2306612 :             IF (k < cyclen) THEN
    1069           0 :                DO l = k + 1, cyclen
    1070           0 :                   rm1(:) = pos(:, p(k), j) - pos(:, p(l), j)
    1071           0 :                   rm2(:) = pos(:, p(k), pk1) - pos(:, p(l), pk1)
    1072           0 :                   ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
    1073           0 :                   rm1(:) = work(:, p(k), j) - work(:, p(l), j)
    1074           0 :                   rm2(:) = work(:, p(k), pk1) - work(:, p(l), pk1)
    1075           0 :                   ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
    1076             :                END DO
    1077             :             END IF
    1078             :          END DO
    1079             :       END DO
    1080      164758 :       j = helium%beads
    1081      164758 :       pk1 = 1
    1082      329516 :       DO k = 1, cyclen
    1083      659032 :          tmp1(:) = pos(:, p(k), j) - pos(:, perm(p(k)), 1)
    1084      164758 :          CALL helium_pbc(helium, tmp1)
    1085      164758 :          ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
    1086      659032 :          tmp1(:) = work(:, p(k), j) - work(:, perm(p(1 + MOD(k, cyclen))), 1)
    1087      164758 :          CALL helium_pbc(helium, tmp1)
    1088      164758 :          ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
    1089     4906761 :          DO l = 1, helium%atoms
    1090     4906761 :             IF (l /= p(k)) THEN
    1091    18308980 :                rm1(:) = pos(:, p(k), j) - pos(:, l, j)
    1092    18308980 :                rm2(:) = pos(:, perm(p(k)), 1) - pos(:, perm(l), 1)
    1093     4577245 :                ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
    1094             :             END IF
    1095             :          END DO
    1096             :          ! counted p[k], p[m] twice. subtract those again
    1097      329516 :          IF (k < cyclen) THEN
    1098           0 :             DO l = k + 1, cyclen
    1099           0 :                rm1(:) = pos(:, p(k), j) - pos(:, p(l), j)
    1100           0 :                rm2(:) = pos(:, perm(p(k)), pk1) - pos(:, perm(p(l)), pk1)
    1101           0 :                ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
    1102             :             END DO
    1103             :          END IF
    1104             :       END DO
    1105      164758 :       IF (cyclen > 1) THEN
    1106             :          !k,c,l
    1107           0 :          c = perm(p(1))
    1108           0 :          DO k = 1, cyclen - 1
    1109           0 :             perm(p(k)) = perm(p(k + 1))
    1110             :          END DO
    1111           0 :          perm(p(cyclen)) = c
    1112             :       END IF
    1113      329516 :       DO k = 1, cyclen
    1114     4906761 :          DO l = 1, helium%atoms
    1115     4906761 :             IF (l /= p(k)) THEN
    1116    18308980 :                rm1(:) = work(:, p(k), j) - work(:, l, j)
    1117    18308980 :                rm2(:) = work(:, perm(p(k)), 1) - work(:, perm(l), 1)
    1118     4577245 :                ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
    1119             :             END IF
    1120             :          END DO
    1121             :          ! counted p[k], p[m] twice. subtract those again
    1122      329516 :          IF (k < cyclen) THEN
    1123           0 :             DO l = k + 1, cyclen
    1124           0 :                rm1(:) = work(:, p(k), j) - work(:, p(l), j)
    1125           0 :                rm2(:) = work(:, perm(p(k)), 1) - work(:, perm(p(l)), 1)
    1126           0 :                ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
    1127             :             END DO
    1128             :          END IF
    1129             :       END DO
    1130      164758 :       DEALLOCATE (work3)
    1131             :       ! ok now accept or reject:
    1132      164758 :       rtmp = helium%rng_stream_uniform%next()
    1133             : !   IF ((dtk+ds-pds<0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
    1134      164758 :       IF (dtk + ds - pds < 0.0_dp) THEN
    1135      100337 :          IF (EXP(dtk + ds - pds) < rtmp) THEN
    1136      320841 :             DO c = ifix, helium%beads
    1137      606033 :                DO k = 1, cyclen
    1138     2281536 :                   work(:, p(k), c) = pos(:, p(k), c)
    1139             :                END DO
    1140             :             END DO
    1141       35649 :             IF (cyclen > 1) THEN
    1142           0 :                c = perm(p(cyclen))
    1143           0 :                DO k = cyclen - 1, 1, -1
    1144           0 :                   perm(p(k + 1)) = perm(p(k))
    1145             :                END DO
    1146           0 :                perm(p(1)) = c
    1147             :             END IF
    1148       35649 :             RETURN
    1149             :          END IF
    1150             :       END IF
    1151             :       ! accepted.
    1152             : 
    1153             :       ! rejection of the whole move if any centroid is farther away
    1154             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    1155             :       ! AND ist not moving towards the center
    1156      129109 :       IF (.NOT. helium%periodic) THEN
    1157       19322 :          IF (helium%solute_present) THEN
    1158       77288 :             new_com(:) = helium%center(:)
    1159             :          ELSE
    1160           0 :             new_com(:) = 0.0_dp
    1161           0 :             DO k = 1, helium%atoms
    1162           0 :                DO l = 1, helium%beads
    1163           0 :                   new_com(:) = new_com(:) + helium%work(:, k, l)
    1164             :                END DO
    1165             :             END DO
    1166           0 :             new_com(:) = new_com(:)/helium%atoms/helium%beads
    1167             :          END IF
    1168             :          ! Cycle through atoms (ignore connectivity)
    1169       19322 :          should_reject = .FALSE.
    1170      114428 :          outer: DO ia = 1, helium%atoms
    1171       95539 :             bis(:) = 0.0_dp
    1172     1624163 :             DO ib = 1, helium%beads
    1173     6210035 :                DO ic = 1, 3
    1174     6114496 :                   bis(ic) = bis(ic) + work(ic, ia, ib) - new_com(ic)
    1175             :                END DO
    1176             :             END DO
    1177      382156 :             bis(:) = bis(:)/helium%beads
    1178      382156 :             rtmp = DOT_PRODUCT(bis, bis)
    1179      114428 :             IF (rtmp .GE. helium%droplet_radius**2) THEN
    1180         433 :                biso(:) = 0.0_dp
    1181        7361 :                DO ib = 1, helium%beads
    1182       28145 :                   DO ic = 1, 3
    1183       27712 :                      biso(ic) = biso(ic) + pos(ic, ia, ib) - helium%center(ic)
    1184             :                   END DO
    1185             :                END DO
    1186        1732 :                biso(:) = biso(:)/helium%beads
    1187        1732 :                rtmpo = DOT_PRODUCT(biso, biso)
    1188             :                ! only reject if it moves away from COM outside the droplet radius
    1189         433 :                IF (rtmpo < rtmp) THEN
    1190             :                   ! found - this one does not fit within R from the center
    1191             :                   should_reject = .TRUE.
    1192             :                   EXIT outer
    1193             :                END IF
    1194             :             END IF
    1195             :          END DO outer
    1196       19322 :          IF (should_reject) THEN
    1197             :             ! restore work and perm, then return
    1198        3897 :             DO c = ifix, helium%beads
    1199        7361 :                DO k = 1, cyclen
    1200       27712 :                   work(:, p(k), c) = pos(:, p(k), c)
    1201             :                END DO
    1202             :             END DO
    1203         433 :             IF (cyclen > 1) THEN
    1204           0 :                c = perm(p(cyclen))
    1205           0 :                DO k = cyclen - 1, 1, -1
    1206           0 :                   perm(p(k + 1)) = perm(p(k))
    1207             :                END DO
    1208           0 :                perm(p(1)) = c
    1209             :             END IF
    1210         433 :             RETURN
    1211             :          END IF
    1212             :       END IF
    1213             :       ! accepted.
    1214             : 
    1215             :       ! copy trial over to the real thing
    1216     1158084 :       DO c = ifix, helium%beads
    1217     2187492 :          DO k = 1, cyclen
    1218     8235264 :             pos(:, p(k), c) = work(:, p(k), c)
    1219             :          END DO
    1220             :       END DO
    1221      257352 :       DO k = 1, cyclen
    1222      257352 :          helium%iperm(perm(p(k))) = p(k)
    1223             :       END DO
    1224      128676 :       helium%num_accepted(level + 2, cyclen) = helium%num_accepted(level + 2, cyclen) + 1
    1225      128676 :       res = .TRUE.
    1226             : 
    1227      128676 :       RETURN
    1228      701096 :    END FUNCTION helium_slice_metro_cyclic
    1229             : 
    1230             : ! ***************************************************************************
    1231             : !> \brief  ...
    1232             : !> \param helium ...
    1233             : !> \param len ...
    1234             : !> \return ...
    1235             : !> \author hforbert
    1236             : ! **************************************************************************************************
    1237      990487 :    FUNCTION helium_select_permutation(helium, len) RESULT(res)
    1238             :       TYPE(helium_solvent_type), POINTER                 :: helium
    1239             :       INTEGER, INTENT(IN)                                :: len
    1240             :       LOGICAL                                            :: res
    1241             : 
    1242             :       INTEGER                                            :: i, j, k, n
    1243      990487 :       INTEGER, DIMENSION(:), POINTER                     :: iperm, p, perm
    1244      990487 :       INTEGER, DIMENSION(:, :), POINTER                  :: nmatrix
    1245             :       REAL(KIND=dp)                                      :: rnd, s1, s2, t
    1246      990487 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ipmatrix, pmatrix, tmatrix
    1247             : 
    1248      990487 :       s1 = 0.0_dp
    1249      990487 :       s2 = 0.0_dp
    1250      990487 :       res = .FALSE.
    1251      990487 :       n = helium%atoms
    1252      990487 :       tmatrix => helium%tmatrix
    1253      990487 :       pmatrix => helium%pmatrix
    1254      990487 :       ipmatrix => helium%ipmatrix
    1255      990487 :       perm => helium%permutation
    1256      990487 :       iperm => helium%iperm
    1257      990487 :       p => helium%ptable
    1258      990487 :       nmatrix => helium%nmatrix
    1259             : 
    1260      990487 :       p(len + 1) = -1
    1261     1980974 :       rnd = helium%rng_stream_uniform%next()
    1262      990487 :       p(1) = INT(n*rnd) + 1
    1263     1028377 :       DO k = 1, len - 1
    1264     1015528 :          t = helium%rng_stream_uniform%next()
    1265             :          ! find the corresponding path to connect to
    1266             :          ! using the precalculated optimal decision tree:
    1267     1015528 :          i = n - 1
    1268             :          DO
    1269     1067754 :             IF (tmatrix(p(k), i) > t) THEN
    1270       51667 :                i = nmatrix(p(k), 2*i - 1)
    1271             :             ELSE
    1272     1016087 :                i = nmatrix(p(k), 2*i)
    1273             :             END IF
    1274     1067754 :             IF (i < 0) EXIT
    1275             :          END DO
    1276     1015528 :          i = -i
    1277             :          ! which particle was it previously connected to?
    1278     1015528 :          p(k + 1) = iperm(i)
    1279             :          ! is it unique? quit if it was already part of the permutation
    1280     1079057 :          DO j = 1, k
    1281     1079057 :             IF (p(j) == p(k + 1)) RETURN
    1282             :          END DO
    1283             :          ! acummulate the needed values for the final
    1284             :          ! accept/reject step:
    1285       37890 :          s1 = s1 + ipmatrix(p(k), i)
    1286       50739 :          s2 = s2 + ipmatrix(p(k), perm(p(k)))
    1287             :       END DO
    1288             :       ! close the permutation loop:
    1289       12849 :       s1 = s1 + ipmatrix(p(len), perm(p(1)))
    1290       12849 :       s2 = s2 + ipmatrix(p(len), perm(p(len)))
    1291             :       ! final accept/reject:
    1292       12849 :       rnd = helium%rng_stream_uniform%next()
    1293       12849 :       t = s1*rnd
    1294       12849 :       IF (t > s2) RETURN
    1295             :       ! ok, we have accepted the permutation
    1296             :       ! calculate the action bias for the subsequent resampling
    1297             :       ! of the paths:
    1298          35 :       s1 = pmatrix(p(len), perm(p(1))) - pmatrix(p(len), perm(p(len)))
    1299          70 :       DO k = 1, len - 1
    1300          70 :          s1 = s1 + pmatrix(p(k), perm(p(k + 1))) - pmatrix(p(k), perm(p(k)))
    1301             :       END DO
    1302          35 :       helium%pweight = s1
    1303          35 :       res = .TRUE.
    1304          35 :       RETURN
    1305      990487 :    END FUNCTION helium_select_permutation
    1306             : 
    1307             : END MODULE helium_sampling

Generated by: LCOV version 1.15