LCOV - code coverage report
Current view: top level - src/motion - helium_sampling.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:dc34ec9) Lines: 560 646 86.7 %
Date: 2023-03-24 20:09:49 Functions: 6 6 100.0 %

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

Generated by: LCOV version 1.15