LCOV - code coverage report
Current view: top level - src/motion - helium_worm.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 1356 1571 86.3 %
Date: 2024-04-20 06:29:22 Functions: 19 19 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 dealing with the canonical worm alogrithm
      10             : !> \author Felix Uhl [fuhl]
      11             : !> \date   2018-10-10
      12             : ! **************************************************************************************************
      13             : MODULE helium_worm
      14             : 
      15             :    USE helium_common,                   ONLY: helium_calc_plength,&
      16             :                                               helium_eval_chain,&
      17             :                                               helium_eval_expansion,&
      18             :                                               helium_pbc,&
      19             :                                               helium_update_coord_system
      20             :    USE helium_interactions,             ONLY: helium_bead_solute_e_f,&
      21             :                                               helium_calc_energy,&
      22             :                                               helium_solute_e_f
      23             :    USE helium_io,                       ONLY: helium_write_line
      24             :    USE helium_types,                    ONLY: helium_solvent_type
      25             :    USE input_constants,                 ONLY: helium_forces_average,&
      26             :                                               helium_forces_last
      27             :    USE kinds,                           ONLY: default_string_length,&
      28             :                                               dp
      29             :    USE mathconstants,                   ONLY: pi
      30             :    USE pint_types,                      ONLY: pint_env_type
      31             : #include "../base/base_uses.f90"
      32             : 
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      37             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_worm'
      38             : 
      39             :    PUBLIC :: helium_sample_worm
      40             : 
      41             : CONTAINS
      42             : 
      43             : ! **************************************************************************************************
      44             : !> \brief Main worm sampling routine
      45             : !> \param helium ...
      46             : !> \param pint_env ...
      47             : !> \author fuhl
      48             : ! **************************************************************************************************
      49          27 :    SUBROUTINE helium_sample_worm(helium, pint_env)
      50             : 
      51             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
      52             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
      53             : 
      54             :       CHARACTER(len=default_string_length)               :: stmp
      55             :       INTEGER :: ac, iatom, ibead, icrawl, iMC, imove, ncentracc, ncentratt, ncloseacc, ncloseatt, &
      56             :          ncrawlbwdacc, ncrawlbwdatt, ncrawlfwdacc, ncrawlfwdatt, ncycle, nMC, nmoveheadacc, &
      57             :          nmoveheadatt, nmovetailacc, nmovetailatt, nopenacc, nopenatt, npswapacc, nstagacc, &
      58             :          nstagatt, nstat, nswapacc, nswapatt, ntot, staging_l
      59             :       REAL(KIND=dp)                                      :: rtmp
      60             : 
      61          27 :       CALL helium_update_coord_system(helium, pint_env)
      62             : 
      63          27 :       ncentratt = 0
      64          27 :       ncentracc = 0
      65          27 :       nstagatt = 0
      66          27 :       nstagacc = 0
      67          27 :       nopenatt = 0
      68          27 :       nopenacc = 0
      69          27 :       ncloseatt = 0
      70          27 :       ncloseacc = 0
      71          27 :       nmoveheadatt = 0
      72          27 :       nmoveheadacc = 0
      73          27 :       nmovetailatt = 0
      74          27 :       nmovetailacc = 0
      75          27 :       ncrawlfwdatt = 0
      76          27 :       ncrawlfwdacc = 0
      77          27 :       ncrawlbwdatt = 0
      78          27 :       ncrawlbwdacc = 0
      79          27 :       nswapatt = 0
      80          27 :       npswapacc = 0
      81          27 :       nswapacc = 0
      82          27 :       nstat = 0
      83          27 :       ntot = 0
      84         108 :       helium%proarea%inst(:) = 0.d0
      85         108 :       helium%prarea2%inst(:) = 0.d0
      86         108 :       helium%wnumber%inst(:) = 0.d0
      87         108 :       helium%wnmber2%inst(:) = 0.d0
      88         108 :       helium%mominer%inst(:) = 0.d0
      89        2844 :       IF (helium%solute_present) helium%force_avrg(:, :) = 0.0d0
      90         297 :       helium%energy_avrg(:) = 0.0d0
      91         831 :       helium%plength_avrg(:) = 0.0d0
      92          27 :       IF (helium%worm_max_open_cycles > 0) THEN
      93           0 :          helium%savepos(:, :, :) = helium%pos(:, :, :)
      94           0 :          helium%saveiperm(:) = helium%iperm(:)
      95           0 :          helium%savepermutation(:) = helium%permutation(:)
      96             :       END IF
      97             : 
      98          27 :       nMC = helium%iter_rot
      99          27 :       ncycle = 0
     100          27 :       IF (helium%worm_allow_open) THEN
     101             :          DO ! Exit criterion at the end of the loop
     102       24622 :             DO iMC = 1, nMC
     103       24510 :                imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
     104       24510 :                IF (helium%worm_is_closed) THEN
     105        6436 :                   IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
     106             :                      ! centroid move
     107         183 :                      iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     108             :                      CALL worm_centroid_move(pint_env, helium, &
     109         183 :                                              iatom, helium%worm_centroid_drmax, ac)
     110         183 :                      ncentratt = ncentratt + 1
     111         183 :                      ncentracc = ncentracc + ac
     112             :                      ! Note: weights for open and centroid move are taken from open sampling
     113             :                      ! staging is adjusted to conserve these weights
     114        6253 :                   ELSE IF ((imove >= helium%worm_centroid_max + 1) .AND. (imove <= helium%worm_open_close_min - 1)) THEN
     115             :                      ! staging move
     116        5461 :                      iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     117        5461 :                      ibead = helium%rng_stream_uniform%next(1, helium%beads)
     118        5461 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     119             :                      CALL worm_staging_move(pint_env, helium, &
     120        5461 :                                             iatom, ibead, staging_l, ac)
     121        5461 :                      nstagatt = nstagatt + 1
     122        5461 :                      nstagacc = nstagacc + ac
     123         792 :                   ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
     124             :                      ! attempt opening of worm
     125         792 :                      iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     126         792 :                      ibead = helium%rng_stream_uniform%next(1, helium%beads)
     127         792 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     128             :                      CALL worm_open_move(pint_env, helium, &
     129         792 :                                          iatom, ibead, staging_l, ac)
     130         792 :                      nopenatt = nopenatt + 1
     131         792 :                      nopenacc = nopenacc + ac
     132             :                   ELSE
     133             :                      ! this must not occur
     134           0 :                      CPABORT("Undefined move selected in helium worm sampling!")
     135             :                   END IF
     136             :                ELSE ! worm is open
     137       18074 :                   IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
     138             :                      ! centroid move
     139         463 :                      iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     140             :                      CALL worm_centroid_move(pint_env, helium, &
     141         463 :                                              iatom, helium%worm_centroid_drmax, ac)
     142         463 :                      ncentratt = ncentratt + 1
     143         463 :                      ncentracc = ncentracc + ac
     144       17611 :                   ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
     145             :                      ! staging move
     146        1047 :                      iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     147        1047 :                      ibead = helium%rng_stream_uniform%next(1, helium%beads)
     148        1047 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     149             :                      CALL worm_staging_move(pint_env, helium, &
     150        1047 :                                             iatom, ibead, staging_l, ac)
     151        1047 :                      nstagatt = nstagatt + 1
     152        1047 :                      nstagacc = nstagacc + ac
     153       16564 :                   ELSE IF ((imove >= helium%worm_fcrawl_min) .AND. (imove <= helium%worm_fcrawl_max)) THEN
     154             :                      ! crawl forward
     155        2820 :                      DO icrawl = 1, helium%worm_repeat_crawl
     156        1880 :                         staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     157             :                         CALL worm_crawl_move_forward(pint_env, helium, &
     158        1880 :                                                      staging_l, ac)
     159        1880 :                         ncrawlfwdatt = ncrawlfwdatt + 1
     160        2820 :                         ncrawlfwdacc = ncrawlfwdacc + ac
     161             :                      END DO
     162       15624 :                   ELSE IF ((imove >= helium%worm_bcrawl_min) .AND. (imove <= helium%worm_bcrawl_max)) THEN
     163             :                      ! crawl backward
     164        3099 :                      DO icrawl = 1, helium%worm_repeat_crawl
     165        2066 :                         staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     166             :                         CALL worm_crawl_move_backward(pint_env, helium, &
     167        2066 :                                                       staging_l, ac)
     168        2066 :                         ncrawlbwdatt = ncrawlbwdatt + 1
     169        3099 :                         ncrawlbwdacc = ncrawlbwdacc + ac
     170             :                      END DO
     171       14591 :                   ELSE IF ((imove >= helium%worm_head_min) .AND. (imove <= helium%worm_head_max)) THEN
     172             :                      ! move head
     173        1033 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     174             :                      CALL worm_head_move(pint_env, helium, &
     175        1033 :                                          staging_l, ac)
     176        1033 :                      nmoveheadatt = nmoveheadatt + 1
     177        1033 :                      nmoveheadacc = nmoveheadacc + ac
     178       13558 :                   ELSE IF ((imove >= helium%worm_tail_min) .AND. (imove <= helium%worm_tail_max)) THEN
     179             :                      ! move tail
     180        1041 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     181             :                      CALL worm_tail_move(pint_env, helium, &
     182        1041 :                                          staging_l, ac)
     183        1041 :                      nmovetailatt = nmovetailatt + 1
     184        1041 :                      nmovetailacc = nmovetailacc + ac
     185       12517 :                   ELSE IF ((imove >= helium%worm_swap_min) .AND. (imove <= helium%worm_swap_max)) THEN
     186       10448 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     187             :                      CALL worm_swap_move(pint_env, helium, &
     188       10448 :                                          helium%atoms, staging_l, ac)
     189       10448 :                      npswapacc = npswapacc + ac
     190       10448 :                      nswapacc = nswapacc + ac
     191       10448 :                      nswapatt = nswapatt + 1
     192        2069 :                   ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
     193             :                      ! attempt closing of worm
     194        2069 :                      staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
     195             :                      CALL worm_close_move(pint_env, helium, &
     196        2069 :                                           staging_l, ac)
     197        2069 :                      ncloseatt = ncloseatt + 1
     198        2069 :                      ncloseacc = ncloseacc + ac
     199             :                   ELSE
     200             :                      ! this must not occur
     201           0 :                      CPABORT("Undefined move selected in helium worm sampling!")
     202             :                   END IF
     203             :                END IF
     204             : 
     205             :                ! Accumulate statistics if we are in the Z-sector:
     206       24510 :                IF (helium%worm_is_closed) THEN
     207        6436 :                   nstat = nstat + 1
     208        6436 :                   IF (helium%solute_present) THEN
     209        6212 :                      IF (helium%get_helium_forces == helium_forces_average) THEN
     210             :                         !TODO needs proper averaging!
     211           0 :                         CALL helium_solute_e_f(pint_env, helium, rtmp)
     212           0 :                         helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
     213             :                      END IF
     214             :                   END IF
     215             :                END IF
     216       24622 :                ntot = ntot + 1
     217             :             END DO ! MC loop
     218             : 
     219         112 :             IF (helium%worm_is_closed) THEN
     220             :                EXIT
     221             :             ELSE
     222          85 :                ncycle = ncycle + 1
     223             :                ! reset position and permutation and start MC cycle again
     224             :                ! if stuck in open position for too long
     225          85 :                IF (helium%worm_max_open_cycles > 0 .AND. ncycle > helium%worm_max_open_cycles) THEN
     226           0 :                   nMC = helium%iter_rot
     227           0 :                   ncycle = 0
     228           0 :                   helium%pos = helium%savepos
     229           0 :                   helium%work = helium%pos
     230           0 :                   helium%permutation = helium%savepermutation
     231           0 :                   helium%iperm = helium%saveiperm
     232           0 :                   CPWARN("Trapped in open state, reset helium to closed state from last MD step.")
     233             :                ELSE
     234          85 :                   nMC = MAX(helium%iter_rot/10, 10)
     235             :                END IF
     236             :             END IF
     237             :          END DO !attempts loop
     238             :       ELSE ! only closed configurations allowed
     239           0 :          DO iMC = 1, nMC
     240           0 :             imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
     241             : 
     242           0 :             IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
     243             :                ! centroid move
     244           0 :                iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     245             :                CALL worm_centroid_move(pint_env, helium, &
     246           0 :                                        iatom, helium%worm_centroid_drmax, ac)
     247           0 :                ncentratt = ncentratt + 1
     248           0 :                ncentracc = ncentracc + ac
     249           0 :             ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
     250             :                ! staging move
     251           0 :                iatom = helium%rng_stream_uniform%next(1, helium%atoms)
     252           0 :                ibead = helium%rng_stream_uniform%next(1, helium%beads)
     253             :                CALL worm_staging_move(pint_env, helium, &
     254           0 :                                       iatom, ibead, helium%worm_staging_l, ac)
     255           0 :                nstagatt = nstagatt + 1
     256           0 :                nstagacc = nstagacc + ac
     257             :             ELSE
     258             :                ! this must not occur
     259           0 :                CPABORT("Undefined move selected in helium worm sampling!")
     260             :             END IF
     261             : 
     262             :             ! Accumulate statistics if we are in closed configurations (which we always are)
     263           0 :             nstat = nstat + 1
     264           0 :             ntot = ntot + 1
     265           0 :             IF (helium%solute_present) THEN
     266           0 :                IF (helium%get_helium_forces == helium_forces_average) THEN
     267             :                   ! TODO: needs proper averaging
     268           0 :                   CALL helium_solute_e_f(pint_env, helium, rtmp)
     269           0 :                   helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
     270             :                END IF
     271             :             END IF
     272             :          END DO ! MC loop
     273             :       END IF
     274             : 
     275             :       ! Save naccepted and ntot
     276             :       helium%num_accepted(1, 1) = ncentracc + nstagacc + nopenacc + ncloseacc + nswapacc + &
     277          27 :                                   nmoveheadacc + nmovetailacc + ncrawlfwdacc + ncrawlbwdacc
     278             :       helium%num_accepted(2, 1) = ncentratt + nstagatt + nopenatt + ncloseatt + nswapatt + &
     279          27 :                                   nmoveheadatt + nmovetailatt + ncrawlfwdatt + ncrawlbwdatt
     280          27 :       helium%worm_nstat = nstat
     281             : 
     282             :       ! Calculate energy and permutation path length
     283             :       ! Multiply times helium%worm_nstat for consistent averaging in helium_sampling
     284          27 :       CALL helium_calc_energy(helium, pint_env)
     285         297 :       helium%energy_avrg(:) = helium%energy_inst(:)*REAL(helium%worm_nstat, dp)
     286             : 
     287          27 :       CALL helium_calc_plength(helium)
     288         831 :       helium%plength_avrg(:) = helium%plength_inst(:)*REAL(helium%worm_nstat, dp)
     289             : 
     290             :       ! Calculate last_force
     291          27 :       IF (helium%solute_present) THEN
     292          17 :          IF (helium%get_helium_forces == helium_forces_last) THEN
     293          17 :             CALL helium_solute_e_f(pint_env, helium, rtmp)
     294        2834 :             helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
     295             :          END IF
     296             :       END IF
     297             : 
     298          27 :       IF (helium%worm_show_statistics) THEN
     299          27 :          WRITE (stmp, *) '--------------------------------------------------'
     300          27 :          CALL helium_write_line(stmp)
     301          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Opening: ', &
     302          27 :             REAL(nopenacc, dp)/REAL(MAX(1, nopenatt), dp), &
     303          54 :             nopenacc, nopenatt
     304          27 :          CALL helium_write_line(stmp)
     305          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Closing: ', &
     306          27 :             REAL(ncloseacc, dp)/REAL(MAX(1, ncloseatt), dp), &
     307          54 :             ncloseacc, ncloseatt
     308          27 :          CALL helium_write_line(stmp)
     309          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Move Head: ', &
     310          27 :             REAL(nmoveheadacc, dp)/REAL(MAX(1, nmoveheadatt), dp), &
     311          54 :             nmoveheadacc, nmoveheadatt
     312          27 :          CALL helium_write_line(stmp)
     313          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Move Tail: ', &
     314          27 :             REAL(nmovetailacc, dp)/REAL(MAX(1, nmovetailatt), dp), &
     315          54 :             nmovetailacc, nmovetailatt
     316          27 :          CALL helium_write_line(stmp)
     317          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Crawl FWD: ', &
     318          27 :             REAL(ncrawlfwdacc, dp)/REAL(MAX(1, ncrawlfwdatt), dp), &
     319          54 :             ncrawlfwdacc, ncrawlfwdatt
     320          27 :          CALL helium_write_line(stmp)
     321          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Crawl BWD: ', &
     322          27 :             REAL(ncrawlbwdacc, dp)/REAL(MAX(1, ncrawlbwdatt), dp), &
     323          54 :             ncrawlbwdacc, ncrawlbwdatt
     324          27 :          CALL helium_write_line(stmp)
     325          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Staging: ', &
     326          27 :             REAL(nstagacc, dp)/REAL(MAX(1, nstagatt), dp), &
     327          54 :             nstagacc, nstagatt
     328          27 :          CALL helium_write_line(stmp)
     329          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Centroid: ', &
     330          27 :             REAL(ncentracc, dp)/REAL(MAX(1, ncentratt), dp), &
     331          54 :             ncentracc, ncentratt
     332          27 :          CALL helium_write_line(stmp)
     333          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Select: ', &
     334          27 :             REAL(npswapacc, dp)/REAL(MAX(1, nswapatt), dp), &
     335          54 :             npswapacc, nswapatt
     336          27 :          CALL helium_write_line(stmp)
     337          27 :          WRITE (stmp, '(A10,F15.5,2I10)') 'Swapping: ', &
     338          27 :             REAL(nswapacc, dp)/REAL(MAX(1, nswapatt), dp), &
     339          54 :             nswapacc, nswapatt
     340          27 :          CALL helium_write_line(stmp)
     341          27 :          WRITE (stmp, *) "Open State Probability:   ", REAL(ntot - nstat, dp)/REAL(MAX(1, ntot), dp), ntot - nstat, ntot
     342          27 :          CALL helium_write_line(stmp)
     343          27 :          WRITE (stmp, *) "Closed State Probability: ", REAL(nstat, dp)/REAL(MAX(1, ntot), dp), nstat, ntot
     344          27 :          CALL helium_write_line(stmp)
     345             :       END IF
     346             : 
     347             :       !CALL center_pos(helium)
     348             : 
     349          27 :    END SUBROUTINE helium_sample_worm
     350             : 
     351             : ! **************************************************************************************************
     352             : !> \brief Centroid Move (accounting for exchanged particles)
     353             : !> \param pint_env ...
     354             : !> \param helium ...
     355             : !> \param iatom ...
     356             : !> \param drmax ...
     357             : !> \param ac ...
     358             : !> \author fuhl
     359             : ! **************************************************************************************************
     360         646 :    SUBROUTINE worm_centroid_move(pint_env, helium, iatom, drmax, ac)
     361             : 
     362             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     363             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     364             :       INTEGER, INTENT(IN)                                :: iatom
     365             :       REAL(KIND=dp), INTENT(IN)                          :: drmax
     366             :       INTEGER, INTENT(OUT)                               :: ac
     367             : 
     368             :       INTEGER                                            :: ia, ib, ic, jatom
     369             :       LOGICAL                                            :: should_reject, worm_in_moved_cycle
     370             :       REAL(KIND=dp)                                      :: rtmp, rtmpo, sdiff, snew, sold
     371             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
     372             : 
     373        2584 :       DO ic = 1, 3
     374        1938 :          rtmp = helium%rng_stream_uniform%next()
     375        2584 :          dr(ic) = (2.0_dp*rtmp - 1.0_dp)*drmax
     376             :       END DO
     377             : 
     378         646 :       IF (helium%worm_is_closed) THEN
     379         183 :          worm_in_moved_cycle = .FALSE.
     380             :          ! Perform move for first atom
     381        3210 :          DO ib = 1, helium%beads
     382       12291 :             helium%work(:, iatom, ib) = helium%work(:, iatom, ib) + dr(:)
     383             :          END DO
     384             :          ! move along permutation cycle
     385         183 :          jatom = helium%permutation(iatom)
     386         183 :          DO WHILE (jatom /= iatom)
     387           0 :             DO ib = 1, helium%beads
     388           0 :                helium%work(:, jatom, ib) = helium%work(:, jatom, ib) + dr(:)
     389             :             END DO
     390             :             ! next atom in chain
     391           0 :             jatom = helium%permutation(jatom)
     392             :          END DO
     393             :       ELSE ! worm is open
     394         463 :          worm_in_moved_cycle = (helium%worm_atom_idx == iatom)
     395             :          ! while moving, check if worm is in moved cycle
     396             :          ! Perform move for first atom
     397        7907 :          DO ib = 1, helium%beads
     398       30239 :             helium%work(:, iatom, ib) = helium%work(:, iatom, ib) + dr(:)
     399             :          END DO
     400             :          ! move along permutation cycle
     401         463 :          jatom = helium%permutation(iatom)
     402         463 :          DO WHILE (jatom /= iatom)
     403           0 :             DO ib = 1, helium%beads
     404           0 :                helium%work(:, jatom, ib) = helium%work(:, jatom, ib) + dr(:)
     405             :             END DO
     406           0 :             worm_in_moved_cycle = worm_in_moved_cycle .OR. (helium%worm_atom_idx == jatom)
     407             :             ! next atom in chain
     408           0 :             jatom = helium%permutation(jatom)
     409             :          END DO
     410             :          ! if atom contains had bead move that as well
     411         463 :          IF (worm_in_moved_cycle) THEN
     412          76 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:) + dr(:)
     413             :          END IF
     414             :       END IF
     415             : 
     416             :       sold = worm_centroid_move_action(helium, helium%pos, iatom, &
     417         646 :                                        helium%worm_xtra_bead, worm_in_moved_cycle)
     418             : 
     419             :       snew = worm_centroid_move_action(helium, helium%work, iatom, &
     420         646 :                                        helium%worm_xtra_bead_work, worm_in_moved_cycle)
     421             : 
     422         646 :       IF (helium%solute_present) THEN
     423             :          sold = sold + worm_centroid_move_inter_action(pint_env, helium, helium%pos, iatom, &
     424         631 :                                                        helium%worm_xtra_bead, worm_in_moved_cycle)
     425             :          snew = snew + worm_centroid_move_inter_action(pint_env, helium, helium%work, iatom, &
     426         631 :                                                        helium%worm_xtra_bead_work, worm_in_moved_cycle)
     427             :       END IF
     428             : 
     429             :       ! Metropolis:
     430         646 :       sdiff = sold - snew
     431         646 :       IF (sdiff < 0) THEN
     432         321 :          should_reject = .FALSE.
     433         321 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
     434             :             should_reject = .TRUE.
     435             :          ELSE
     436         321 :             rtmp = helium%rng_stream_uniform%next()
     437         321 :             IF (EXP(sdiff) < rtmp) THEN
     438             :                should_reject = .TRUE.
     439             :             END IF
     440             :          END IF
     441             :          IF (should_reject) THEN
     442             :             ! rejected !
     443             :             ! write back only changed atoms
     444         206 :             DO ib = 1, helium%beads
     445         794 :                helium%work(:, iatom, ib) = helium%pos(:, iatom, ib)
     446             :             END DO
     447             :             ! move along permutation cycle
     448          10 :             jatom = helium%permutation(iatom)
     449          10 :             DO WHILE (jatom /= iatom)
     450           0 :                DO ib = 1, helium%beads
     451           0 :                   helium%work(:, jatom, ib) = helium%pos(:, jatom, ib)
     452             :                END DO
     453             :                ! next atom in chain
     454           0 :                jatom = helium%permutation(jatom)
     455             :             END DO
     456          40 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
     457          10 :             ac = 0
     458          12 :             RETURN
     459             :          END IF
     460             :       END IF
     461             : 
     462             :       ! for now accepted
     463             :       ! rejection of the whole move if any centroid is farther away
     464             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
     465             :       ! AND ist not moving towards the center
     466         636 :       IF (.NOT. helium%periodic) THEN
     467         625 :          IF (helium%solute_present) THEN
     468        2500 :             new_com(:) = helium%center(:)
     469         625 :             old_com(:) = new_com(:)
     470             :          ELSE
     471           0 :             new_com(:) = 0.0_dp
     472           0 :             DO ia = 1, helium%atoms
     473           0 :                DO ib = 1, helium%beads
     474           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
     475             :                END DO
     476             :             END DO
     477           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
     478             :             ! also compute the old center of mass (optimization potential)
     479           0 :             old_com(:) = 0.0_dp
     480           0 :             DO ia = 1, helium%atoms
     481           0 :                DO ib = 1, helium%beads
     482           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
     483             :                END DO
     484             :             END DO
     485           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
     486             :          END IF
     487         625 :          should_reject = .FALSE.
     488       18905 :          atom: DO ia = 1, helium%atoms
     489       18282 :             dr(:) = 0.0_dp
     490      310794 :             DO ib = 1, helium%beads
     491     1188330 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
     492             :             END DO
     493       73128 :             dr(:) = dr(:)/REAL(helium%beads, dp)
     494       73128 :             rtmp = DOT_PRODUCT(dr, dr)
     495       18905 :             IF (rtmp >= helium%droplet_radius**2) THEN
     496           2 :                dro(:) = 0.0_dp
     497          34 :                DO ib = 1, helium%beads
     498         130 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
     499             :                END DO
     500           8 :                dro(:) = dro(:)/REAL(helium%beads, dp)
     501           8 :                rtmpo = DOT_PRODUCT(dro, dro)
     502             :                ! only reject if it moves away from COM outside the droplet radius
     503           2 :                IF (rtmpo < rtmp) THEN
     504             :                   ! found - this one does not fit within R from the center
     505             :                   should_reject = .TRUE.
     506             :                   EXIT atom
     507             :                END IF
     508             :             END IF
     509             :          END DO atom
     510         625 :          IF (should_reject) THEN
     511             :             ! restore original coordinates
     512             :             ! write back only changed atoms
     513          34 :             DO ib = 1, helium%beads
     514         130 :                helium%work(:, iatom, ib) = helium%pos(:, iatom, ib)
     515             :             END DO
     516             :             ! move along permutation cycle
     517           2 :             jatom = helium%permutation(iatom)
     518           2 :             DO WHILE (jatom /= iatom)
     519           0 :                DO ib = 1, helium%beads
     520           0 :                   helium%work(:, jatom, ib) = helium%pos(:, jatom, ib)
     521             :                END DO
     522             :                ! next atom in chain
     523           0 :                jatom = helium%permutation(jatom)
     524             :             END DO
     525           8 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
     526           2 :             ac = 0
     527           2 :             RETURN
     528             :          END IF
     529             :       END IF
     530             : 
     531             :       ! finally accepted
     532         634 :       ac = 1
     533             :       ! write changed coordinates to position array
     534       10877 :       DO ib = 1, helium%beads
     535       41606 :          helium%pos(:, iatom, ib) = helium%work(:, iatom, ib)
     536             :       END DO
     537             :       ! move along permutation cycle
     538         634 :       jatom = helium%permutation(iatom)
     539         634 :       DO WHILE (jatom /= iatom)
     540           0 :          DO ib = 1, helium%beads
     541           0 :             helium%pos(:, jatom, ib) = helium%work(:, jatom, ib)
     542             :          END DO
     543             :          ! next atom in chain
     544           0 :          jatom = helium%permutation(jatom)
     545             :       END DO
     546        2536 :       helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
     547             : 
     548             :    END SUBROUTINE worm_centroid_move
     549             : 
     550             : ! **************************************************************************************************
     551             : !> \brief Action for centroid move
     552             : !> \param helium ...
     553             : !> \param pos ...
     554             : !> \param iatom ...
     555             : !> \param xtrapos ...
     556             : !> \param winc ...
     557             : !> \return ...
     558             : !> \author fuhl
     559             : ! **************************************************************************************************
     560        1292 :    REAL(KIND=dp) FUNCTION worm_centroid_move_action(helium, pos, iatom, xtrapos, winc) &
     561             :       RESULT(partaction)
     562             : 
     563             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     564             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
     565             :          POINTER                                         :: pos
     566             :       INTEGER, INTENT(IN)                                :: iatom
     567             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: xtrapos
     568             :       LOGICAL, INTENT(IN)                                :: winc
     569             : 
     570             :       INTEGER                                            :: ia, ib, jatom, katom, opatom, patom, &
     571             :                                                             wbead
     572             :       LOGICAL                                            :: incycle
     573        1292 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work2, work3
     574        1292 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     575             :       REAL(KIND=dp), DIMENSION(3)                        :: r, rp
     576             : 
     577        9044 :       ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
     578             : 
     579             :       ! compute action difference
     580             :       ! action before the move from pos array
     581        1292 :       partaction = 0.0_dp
     582             : 
     583             :       ! pair of first atom with non-in-cycle atoms
     584        1292 :       jatom = iatom
     585       39204 :       DO ia = 1, helium%atoms
     586       37912 :          IF (ia == jatom) CYCLE
     587       36620 :          incycle = .FALSE.
     588       36620 :          katom = helium%permutation(jatom)
     589       36620 :          DO WHILE (katom /= jatom)
     590           0 :             IF (katom == ia) THEN
     591             :                incycle = .TRUE.
     592             :                EXIT
     593             :             END IF
     594           0 :             katom = helium%permutation(katom)
     595             :          END DO
     596       36620 :          IF (incycle) CYCLE
     597             :          ! if not in cycle, compute pair action
     598      630910 :          DO ib = 1, helium%beads
     599     2413780 :             work(:, ib) = pos(:, ia, ib) - pos(:, jatom, ib)
     600             :          END DO
     601             :          work(:, helium%beads + 1) = pos(:, helium%permutation(ia), 1) - &
     602      146480 :                                      pos(:, helium%permutation(jatom), 1)
     603       39204 :          partaction = partaction + helium_eval_chain(helium, work, helium%beads + 1, work2, work3)
     604             :       END DO
     605             :       ! all other cycle atoms with non-in-cycle atoms
     606        1292 :       jatom = helium%permutation(iatom)
     607        1292 :       DO WHILE (jatom /= iatom)
     608           0 :          DO ia = 1, helium%atoms
     609           0 :             IF (ia == jatom) CYCLE
     610           0 :             incycle = .FALSE.
     611           0 :             katom = helium%permutation(jatom)
     612           0 :             DO WHILE (katom /= jatom)
     613           0 :                IF (katom == ia) THEN
     614             :                   incycle = .TRUE.
     615             :                   EXIT
     616             :                END IF
     617           0 :                katom = helium%permutation(katom)
     618             :             END DO
     619           0 :             IF (incycle) CYCLE
     620             :             ! if not in cycle, compute pair action
     621           0 :             DO ib = 1, helium%beads
     622           0 :                work(:, ib) = pos(:, ia, ib) - pos(:, jatom, ib)
     623             :             END DO
     624             :             work(:, helium%beads + 1) = pos(:, helium%permutation(ia), 1) - &
     625           0 :                                         pos(:, helium%permutation(jatom), 1)
     626           0 :             partaction = partaction + helium_eval_chain(helium, work, helium%beads + 1, work2, work3)
     627             :          END DO
     628           0 :          jatom = helium%permutation(jatom)
     629             :       END DO
     630             :       ! correct pair action for open worm configurations
     631        1292 :       IF (.NOT. helium%worm_is_closed) THEN
     632         926 :          wbead = helium%worm_bead_idx
     633         926 :          IF (winc) THEN
     634          38 :             IF (helium%worm_bead_idx == 1) THEN
     635             :                ! patom is the atom in front of the lone head bead
     636           8 :                patom = helium%iperm(helium%worm_atom_idx)
     637             :                ! go through all atoms, not in the cycle, and correct pair action
     638         264 :                DO ia = 1, helium%atoms
     639         256 :                   IF (ia == helium%worm_atom_idx) CYCLE
     640         248 :                   incycle = .FALSE.
     641         248 :                   katom = helium%permutation(helium%worm_atom_idx)
     642         248 :                   DO WHILE (katom /= helium%worm_atom_idx)
     643           0 :                      IF (katom == ia) THEN
     644             :                         incycle = .TRUE.
     645             :                         EXIT
     646             :                      END IF
     647           0 :                      katom = helium%permutation(katom)
     648             :                   END DO
     649         248 :                   IF (incycle) CYCLE
     650             :                   ! find other previous atom
     651             :                   ! opatom is the atom in front of the first bead of the second atom
     652         248 :                   opatom = helium%iperm(ia)
     653             :                   ! if not in cycle, compute pair action
     654             :                   ! subtract pair action for closed link
     655         992 :                   r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, ia, 1)
     656         992 :                   rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
     657         248 :                   partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
     658             :                   ! and add corrected extra link
     659             :                   ! rp stays the same
     660         992 :                   r(:) = xtrapos(:) - pos(:, ia, 1)
     661         264 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
     662             :                END DO
     663             :             ELSE ! bead index /= 1
     664             :                ! atom index is constant
     665             :                ! go through all atoms, not in the cycle, and correct pair action
     666         870 :                DO ia = 1, helium%atoms
     667         840 :                   IF (ia == helium%worm_atom_idx) CYCLE
     668         810 :                   incycle = .FALSE.
     669         810 :                   katom = helium%permutation(helium%worm_atom_idx)
     670         810 :                   DO WHILE (katom /= helium%worm_atom_idx)
     671           0 :                      IF (katom == ia) THEN
     672             :                         incycle = .TRUE.
     673             :                         EXIT
     674             :                      END IF
     675           0 :                      katom = helium%permutation(katom)
     676             :                   END DO
     677         810 :                   IF (incycle) CYCLE
     678             :                   ! if not in cycle, compute pair action
     679             :                   ! subtract pair action for closed link
     680        3240 :                   r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, ia, wbead)
     681        3240 :                   rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, ia, wbead - 1)
     682         810 :                   partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
     683             :                   ! and add corrected extra link
     684             :                   ! rp stays the same
     685        3240 :                   r(:) = xtrapos(:) - pos(:, ia, wbead)
     686         870 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
     687             :                END DO
     688             :             END IF
     689             :          ELSE ! worm is not the moved cycle
     690         888 :             IF (helium%worm_bead_idx == 1) THEN
     691             :                ! patom is the atom in front of the lone head bead
     692          60 :                patom = helium%iperm(helium%worm_atom_idx)
     693          60 :                opatom = helium%iperm(iatom)
     694             :                !correct action contribution for first atom in moved cycle
     695         240 :                r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, iatom, 1)
     696         240 :                rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
     697          60 :                partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
     698             :                ! and add corrected extra link
     699             :                ! rp stays the same
     700         240 :                r(:) = xtrapos(:) - pos(:, iatom, 1)
     701          60 :                partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
     702             :                ! go through all other atoms, not in the exchange cycle, and correct pair action
     703          60 :                ia = helium%permutation(iatom)
     704          60 :                DO WHILE (ia /= iatom)
     705           0 :                   opatom = helium%iperm(ia)
     706             :                   ! subtract pair action for closed link
     707           0 :                   r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, ia, 1)
     708           0 :                   rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
     709           0 :                   partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
     710             :                   ! and add corrected extra link
     711             :                   ! rp stays the same
     712           0 :                   r(:) = xtrapos(:) - pos(:, ia, 1)
     713           0 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
     714          60 :                   ia = helium%permutation(ia)
     715             :                END DO
     716             :             ELSE ! bead index /= 1
     717             :                ! patom is the atom in front of the lone head bead
     718             :                !correct action contribution for first atom in moved cycle
     719        3312 :                r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, iatom, wbead)
     720        3312 :                rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, iatom, wbead - 1)
     721         828 :                partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
     722             :                ! and add corrected extra link
     723             :                ! rp stays the same
     724        3312 :                r(:) = xtrapos(:) - pos(:, iatom, wbead)
     725         828 :                partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
     726             :                ! go through all other atoms, not in the exchange cycle, and correct pair action
     727         828 :                ia = helium%permutation(iatom)
     728         828 :                DO WHILE (ia /= iatom)
     729             :                   ! subtract pair action for closed link
     730           0 :                   r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, ia, wbead)
     731           0 :                   rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, ia, wbead - 1)
     732           0 :                   partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
     733             :                   ! and add corrected extra link
     734             :                   ! rp stays the same
     735           0 :                   r(:) = xtrapos(:) - pos(:, ia, wbead)
     736           0 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
     737         828 :                   ia = helium%permutation(ia)
     738             :                END DO
     739             :             END IF
     740             :          END IF
     741             :       END IF
     742        1292 :       DEALLOCATE (work, work2, work3)
     743             : 
     744        1292 :    END FUNCTION worm_centroid_move_action
     745             : 
     746             : ! **************************************************************************************************
     747             : !> \brief Centroid interaction
     748             : !> \param pint_env ...
     749             : !> \param helium ...
     750             : !> \param pos ...
     751             : !> \param iatom ...
     752             : !> \param xtrapos ...
     753             : !> \param winc ...
     754             : !> \return ...
     755             : !> \author fuhl
     756             : ! **************************************************************************************************
     757        1262 :    REAL(KIND=dp) FUNCTION worm_centroid_move_inter_action(pint_env, helium, pos, iatom, xtrapos, winc) &
     758             :       RESULT(partaction)
     759             : 
     760             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     761             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     762             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
     763             :          POINTER                                         :: pos
     764             :       INTEGER, INTENT(IN)                                :: iatom
     765             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: xtrapos
     766             :       LOGICAL, INTENT(IN)                                :: winc
     767             : 
     768             :       INTEGER                                            :: jatom, jbead
     769             :       REAL(KIND=dp)                                      :: energy
     770             : 
     771        1262 :       partaction = 0.0_dp
     772             :       ! spcial treatment if worm is in moved cycle
     773        1262 :       IF (winc) THEN
     774             :          ! first atom by hand
     775          38 :          jatom = iatom
     776             :          ! if it is worm atom it gets special treatment
     777          38 :          IF (jatom == helium%worm_atom_idx) THEN
     778             :             ! up to worm intersection
     779         260 :             DO jbead = 1, helium%worm_bead_idx - 1
     780             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     781         222 :                                            jatom, jbead, pos(:, jatom, jbead), energy=energy)
     782         260 :                partaction = partaction + energy
     783             :             END DO
     784             :             ! head and tail each with 1/2 weight
     785          38 :             jbead = helium%worm_bead_idx
     786             :             ! tail
     787             :             CALL helium_bead_solute_e_f(pint_env, helium, &
     788          38 :                                         jatom, jbead, pos(:, jatom, jbead), energy=energy)
     789          38 :             partaction = partaction + 0.5_dp*energy
     790             :             ! head
     791             :             CALL helium_bead_solute_e_f(pint_env, helium, &
     792          38 :                                         jatom, jbead, xtrapos, energy=energy)
     793          38 :             partaction = partaction + 0.5_dp*energy
     794             :             ! rest of ring polymer
     795         386 :             DO jbead = helium%worm_bead_idx + 1, helium%beads
     796             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     797         348 :                                            jatom, jbead, pos(:, jatom, jbead), energy=energy)
     798         386 :                partaction = partaction + energy
     799             :             END DO
     800             :          ELSE
     801           0 :             DO jbead = 1, helium%beads
     802             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     803           0 :                                            jatom, jbead, pos(:, jatom, jbead), energy=energy)
     804           0 :                partaction = partaction + energy
     805             :             END DO
     806             :          END IF
     807          38 :          jatom = helium%permutation(iatom)
     808          38 :          DO WHILE (jatom /= iatom)
     809           0 :             IF (jatom == helium%worm_atom_idx) THEN
     810             :                ! up to worm intersection
     811           0 :                DO jbead = 1, helium%worm_bead_idx - 1
     812             :                   CALL helium_bead_solute_e_f(pint_env, helium, &
     813           0 :                                               jatom, jbead, pos(:, jatom, jbead), energy=energy)
     814           0 :                   partaction = partaction + energy
     815             :                END DO
     816             :                ! head and tail each with 1/2 weight
     817           0 :                jbead = helium%worm_bead_idx
     818             :                ! tail
     819             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     820           0 :                                            jatom, jbead, pos(:, jatom, jbead), energy=energy)
     821           0 :                partaction = partaction + 0.5_dp*energy
     822             :                ! head
     823             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     824           0 :                                            jatom, jbead, xtrapos, energy=energy)
     825           0 :                partaction = partaction + 0.5_dp*energy
     826             :                ! rest of ring polymer
     827           0 :                DO jbead = helium%worm_bead_idx + 1, helium%beads
     828             :                   CALL helium_bead_solute_e_f(pint_env, helium, &
     829           0 :                                               jatom, jbead, pos(:, jatom, jbead), energy=energy)
     830           0 :                   partaction = partaction + energy
     831             :                END DO
     832             :             ELSE
     833           0 :                DO jbead = 1, helium%beads
     834             :                   CALL helium_bead_solute_e_f(pint_env, helium, &
     835           0 :                                               jatom, jbead, pos(:, jatom, jbead), energy=energy)
     836           0 :                   partaction = partaction + energy
     837             :                END DO
     838             :             END IF
     839           0 :             jatom = helium%permutation(jatom)
     840             :          END DO
     841             :       ELSE ! worm not in moved cycle
     842             :          ! first atom by hand
     843        1224 :          jatom = iatom
     844             :          ! if it is worm atom it gets special treatment
     845       20808 :          DO jbead = 1, helium%beads
     846             :             CALL helium_bead_solute_e_f(pint_env, helium, &
     847       19584 :                                         jatom, jbead, pos(:, jatom, jbead), energy=energy)
     848       20808 :             partaction = partaction + energy
     849             :          END DO
     850        1224 :          jatom = helium%permutation(iatom)
     851        1224 :          DO WHILE (jatom /= iatom)
     852           0 :             DO jbead = 1, helium%beads
     853             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     854           0 :                                            jatom, jbead, pos(:, jatom, jbead), energy=energy)
     855           0 :                partaction = partaction + energy
     856             :             END DO
     857        1224 :             jatom = helium%permutation(jatom)
     858             :          END DO
     859             :       END IF
     860             : 
     861        1262 :       partaction = partaction*helium%tau
     862             : 
     863        1262 :    END FUNCTION worm_centroid_move_inter_action
     864             : 
     865             : ! **************************************************************************************************
     866             : !> \brief General path construct based on staging moves
     867             : !> \param helium ...
     868             : !> \param ri ...
     869             : !> \param rj ...
     870             : !> \param l ...
     871             : !> \param new_path ...
     872             : !> \author fuhl
     873             : ! **************************************************************************************************
     874       19011 :    SUBROUTINE path_construct(helium, ri, rj, l, new_path)
     875             : 
     876             :       !constructs a path of length l between the positions ri and rj
     877             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     878             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: ri, rj
     879             :       INTEGER, INTENT(IN)                                :: l
     880             :       REAL(KIND=dp), DIMENSION(3, l), INTENT(OUT)        :: new_path
     881             : 
     882             :       INTEGER                                            :: idim, istage
     883             :       REAL(KIND=dp)                                      :: imass, invstagemass, rk, weight
     884             :       REAL(KIND=dp), DIMENSION(3)                        :: re, rs
     885             : 
     886       19011 :       imass = 1.0_dp/helium%he_mass_au
     887             :       ! dealing with periodicity
     888       19011 :       rs(:) = ri(:)
     889       76044 :       re(:) = rj(:) - rs(:)
     890       19011 :       CALL helium_pbc(helium, re)
     891       76044 :       re(:) = re(:) + rs(:)
     892             : 
     893             :       ! first construction by hand
     894             :       ! reusable weight factor 1/(l+1)
     895       19011 :       rk = REAL(l, dp)
     896       19011 :       weight = 1.0_dp/(rk + 1.0_dp)
     897             :       ! staging mass needed for modified variance
     898       19011 :       invstagemass = rk*weight*imass
     899             :       ! proposing new positions
     900       76044 :       DO idim = 1, 3
     901       76044 :          new_path(idim, 1) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
     902             :       END DO
     903       76044 :       new_path(:, 1) = new_path(:, 1) + weight*(re(:) + rk*rs(:))
     904             : 
     905       64566 :       DO istage = 2, l
     906             :          ! reusable weight factor 1/(k+1)
     907       45555 :          rk = REAL(l - istage + 1, dp)
     908       45555 :          weight = 1.0_dp/(rk + 1.0_dp)
     909             :          ! staging mass needed for modified variance
     910       45555 :          invstagemass = rk*weight*imass
     911             :          ! proposing new positions
     912      182220 :          DO idim = 1, 3
     913      182220 :             new_path(idim, istage) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
     914             :          END DO
     915      201231 :          new_path(:, istage) = new_path(:, istage) + weight*(rk*new_path(:, istage - 1) + re(:))
     916             :       END DO
     917             : 
     918       19011 :    END SUBROUTINE path_construct
     919             : 
     920             : ! **************************************************************************************************
     921             : !> \brief Staging move
     922             : !> \param pint_env ...
     923             : !> \param helium ...
     924             : !> \param startatom ...
     925             : !> \param startbead ...
     926             : !> \param l ...
     927             : !> \param ac ...
     928             : !> \author fuhl
     929             : ! **************************************************************************************************
     930        6508 :    SUBROUTINE worm_staging_move(pint_env, helium, startatom, startbead, l, ac)
     931             : 
     932             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     933             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     934             :       INTEGER, INTENT(IN)                                :: startatom, startbead, l
     935             :       INTEGER, INTENT(OUT)                               :: ac
     936             : 
     937             :       INTEGER                                            :: endatom, endbead, ia, ib, ibead, jbead
     938             :       LOGICAL                                            :: should_reject, worm_overlap
     939             :       REAL(KIND=dp)                                      :: rtmp, rtmpo, sdiff, snew, sold
     940             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
     941        6508 :       REAL(KIND=dp), DIMENSION(3, l)                     :: newsection
     942             : 
     943        6508 :       ac = 0
     944        6508 :       endbead = startbead + l + 1
     945             :       ! Check if the imaginary time section belongs to two atoms
     946        6508 :       IF (endbead > helium%beads) THEN
     947        1833 :          endatom = helium%permutation(startatom)
     948        1833 :          endbead = endbead - helium%beads
     949             :       ELSE
     950        4675 :          endatom = startatom
     951             :       END IF
     952             : 
     953             :       ! check if the imaginary time section overlaps with the worm opening
     954        6508 :       IF (helium%worm_is_closed) THEN
     955             :          worm_overlap = .FALSE.
     956             :       ELSE
     957             :          ! if it does give it special treatment during action evaluation
     958             :          worm_overlap = ((startbead < endbead) .AND. &
     959             :                          (helium%worm_bead_idx > startbead) .AND. &
     960             :                          (helium%worm_bead_idx <= endbead)) &
     961             :                         .OR. &
     962             :                         ((startbead > endbead) .AND. &
     963             :                          ((helium%worm_bead_idx > startbead) .OR. &
     964        1047 :                           (helium%worm_bead_idx <= endbead)))
     965             :          IF (worm_overlap) THEN
     966             :             ! if in addition the worm end is IN the reconstructed section reject immediately
     967         317 :             IF (((helium%worm_atom_idx == startatom) .AND. (helium%worm_bead_idx >= startbead)) .OR. &
     968             :                 ((helium%worm_atom_idx == endatom) .AND. (helium%worm_bead_idx <= endbead))) THEN
     969             :                ac = 0
     970             :                RETURN
     971             :             END IF
     972             :          END IF
     973             :       END IF
     974             :       ! action before the move
     975        1033 :       IF (worm_overlap) THEN
     976             :          sold = worm_path_action_worm_corrected(helium, helium%pos, &
     977             :                                                 startatom, startbead, endatom, endbead, &
     978         303 :                                                 helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
     979             :       ELSE
     980             :          sold = worm_path_action(helium, helium%pos, &
     981        6191 :                                  startatom, startbead, endatom, endbead)
     982             :       END IF
     983             : 
     984        6494 :       IF (helium%solute_present) THEN
     985             :          ! no special head treatment needed, because a swap can't go over
     986             :          ! the worm gap and due to primitive coupling no cross bead terms are considered
     987             :          sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
     988        6291 :                                               startatom, startbead, endatom, endbead)
     989             :       END IF
     990             : 
     991             :       ! construct a new path connecting the start and endbead
     992             :       CALL path_construct(helium, &
     993             :                           helium%pos(:, startatom, startbead), &
     994             :                           helium%pos(:, endatom, endbead), l, &
     995        6494 :                           newsection)
     996             : 
     997             :       ! write new path segment to work array
     998             :       ! first the part that is guaranteed to fit on the coorinates of startatom
     999        6494 :       jbead = 1
    1000       25813 :       DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
    1001       77276 :          helium%work(:, startatom, ibead) = newsection(:, jbead)
    1002       25813 :          jbead = jbead + 1
    1003             :       END DO
    1004             :       ! transfer the rest of the beads to coordinates of endatom if necessary
    1005        6494 :       IF (helium%beads < startbead + l) THEN
    1006        4918 :          DO ibead = 1, endbead - 1
    1007       13956 :             helium%work(:, endatom, ibead) = newsection(:, jbead)
    1008        4918 :             jbead = jbead + 1
    1009             :          END DO
    1010             :       END IF
    1011             : 
    1012             :       ! action after the move
    1013        6494 :       IF (worm_overlap) THEN
    1014             :          snew = worm_path_action_worm_corrected(helium, helium%work, &
    1015             :                                                 startatom, startbead, endatom, endbead, &
    1016         303 :                                                 helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    1017             :       ELSE
    1018             :          snew = worm_path_action(helium, helium%work, &
    1019        6191 :                                  startatom, startbead, endatom, endbead)
    1020             :       END IF
    1021             : 
    1022        6494 :       IF (helium%solute_present) THEN
    1023             :          ! no special head treatment needed, because a swap can't go over
    1024             :          ! the worm gap and due to primitive coupling no cross bead terms are considered
    1025             :          snew = snew + worm_path_inter_action(pint_env, helium, helium%work, &
    1026        6291 :                                               startatom, startbead, endatom, endbead)
    1027             :       END IF
    1028             : 
    1029             :       ! Metropolis:
    1030        6494 :       sdiff = sold - snew
    1031        6494 :       IF (sdiff < 0) THEN
    1032        3303 :          should_reject = .FALSE.
    1033        3303 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    1034             :             should_reject = .TRUE.
    1035             :          ELSE
    1036        3302 :             rtmp = helium%rng_stream_uniform%next()
    1037        3302 :             IF (EXP(sdiff) < rtmp) THEN
    1038             :                should_reject = .TRUE.
    1039             :             END IF
    1040             :          END IF
    1041             :          IF (should_reject) THEN
    1042             :             ! rejected !
    1043             :             ! write back only changed atoms
    1044          98 :             jbead = 1
    1045         446 :             DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
    1046        1392 :                helium%work(:, startatom, ibead) = helium%pos(:, startatom, ibead)
    1047         446 :                jbead = jbead + 1
    1048             :             END DO
    1049             :             ! transfer the rest of the beads to coordinates of endatom if necessary
    1050          98 :             IF (helium%beads < startbead + l) THEN
    1051          54 :                DO ibead = 1, endbead - 1
    1052         148 :                   helium%work(:, endatom, ibead) = helium%pos(:, endatom, ibead)
    1053          54 :                   jbead = jbead + 1
    1054             :                END DO
    1055             :             END IF
    1056          98 :             ac = 0
    1057          98 :             RETURN
    1058             :          END IF
    1059             :       END IF
    1060             : 
    1061             :       ! for now accepted
    1062             :       ! rejection of the whole move if any centroid is farther away
    1063             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    1064             :       ! AND ist not moving towards the center
    1065        6396 :       IF (.NOT. helium%periodic) THEN
    1066        6231 :          IF (helium%solute_present) THEN
    1067       24924 :             new_com(:) = helium%center(:)
    1068        6231 :             old_com(:) = new_com(:)
    1069             :          ELSE
    1070           0 :             new_com(:) = 0.0_dp
    1071           0 :             DO ia = 1, helium%atoms
    1072           0 :                DO ib = 1, helium%beads
    1073           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    1074             :                END DO
    1075             :             END DO
    1076           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1077             :             ! also compute the old center of mass (optimization potential)
    1078           0 :             old_com(:) = 0.0_dp
    1079           0 :             DO ia = 1, helium%atoms
    1080           0 :                DO ib = 1, helium%beads
    1081           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    1082             :                END DO
    1083             :             END DO
    1084           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1085             :          END IF
    1086        6231 :          should_reject = .FALSE.
    1087      182841 :          atom: DO ia = 1, helium%atoms
    1088      176640 :             dr(:) = 0.0_dp
    1089     3002880 :             DO ib = 1, helium%beads
    1090    11481600 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    1091             :             END DO
    1092      706560 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    1093      706560 :             rtmp = DOT_PRODUCT(dr, dr)
    1094      182841 :             IF (rtmp >= helium%droplet_radius**2) THEN
    1095          30 :                dro(:) = 0.0_dp
    1096         510 :                DO ib = 1, helium%beads
    1097        1950 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    1098             :                END DO
    1099         120 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    1100         120 :                rtmpo = DOT_PRODUCT(dro, dro)
    1101             :                ! only reject if it moves away from COM outside the droplet radius
    1102          30 :                IF (rtmpo < rtmp) THEN
    1103             :                   ! found - this one does not fit within R from the center
    1104             :                   should_reject = .TRUE.
    1105             :                   EXIT atom
    1106             :                END IF
    1107             :             END IF
    1108             :          END DO atom
    1109        6231 :          IF (should_reject) THEN
    1110             :             ! restore original coordinates
    1111             :             ! write back only changed atoms
    1112          30 :             jbead = 1
    1113         107 :             DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
    1114         308 :                helium%work(:, startatom, ibead) = helium%pos(:, startatom, ibead)
    1115         107 :                jbead = jbead + 1
    1116             :             END DO
    1117             :             ! transfer the rest of the beads to coordinates of endatom if necessary
    1118          30 :             IF (helium%beads < startbead + l) THEN
    1119          42 :                DO ibead = 1, endbead - 1
    1120         124 :                   helium%work(:, endatom, ibead) = helium%pos(:, endatom, ibead)
    1121          42 :                   jbead = jbead + 1
    1122             :                END DO
    1123             :             END IF
    1124          30 :             ac = 0
    1125          30 :             RETURN
    1126             :          END IF
    1127             :       END IF
    1128             : 
    1129             :       ! finally accepted
    1130        6366 :       ac = 1
    1131             :       ! write changed coordinates to position array
    1132        6366 :       jbead = 1
    1133       25260 :       DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
    1134       75576 :          helium%pos(:, startatom, ibead) = helium%work(:, startatom, ibead)
    1135       25260 :          jbead = jbead + 1
    1136             :       END DO
    1137             :       ! transfer the rest of the beads to coordinates of endatom if necessary
    1138        6366 :       IF (helium%beads < startbead + l) THEN
    1139        4822 :          DO ibead = 1, endbead - 1
    1140       13684 :             helium%pos(:, endatom, ibead) = helium%work(:, endatom, ibead)
    1141        4822 :             jbead = jbead + 1
    1142             :          END DO
    1143             :       END IF
    1144             : 
    1145             :    END SUBROUTINE worm_staging_move
    1146             : 
    1147             : ! **************************************************************************************************
    1148             : !> \brief Open move to off-diagonal elements of the density matrix, allows to sample permutations
    1149             : !> \param pint_env ...
    1150             : !> \param helium ...
    1151             : !> \param endatom ...
    1152             : !> \param endbead ...
    1153             : !> \param l ...
    1154             : !> \param ac ...
    1155             : !> \author fuhl
    1156             : ! **************************************************************************************************
    1157         792 :    SUBROUTINE worm_open_move(pint_env, helium, endatom, endbead, l, ac)
    1158             : 
    1159             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    1160             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1161             :       INTEGER, INTENT(IN)                                :: endatom, endbead, l
    1162             :       INTEGER, INTENT(OUT)                               :: ac
    1163             : 
    1164             :       INTEGER                                            :: ia, ib, idim, kbead, startatom, startbead
    1165             :       LOGICAL                                            :: should_reject
    1166             :       REAL(KIND=dp)                                      :: distsq, mass, rtmp, rtmpo, sdiff, snew, &
    1167             :                                                             sold, xr
    1168             :       REAL(KIND=dp), DIMENSION(3)                        :: distvec, dr, dro, new_com, old_com
    1169             : 
    1170         792 :       mass = helium%he_mass_au
    1171             : 
    1172             :       ! get index of the atom and bead, where the resampling of the head begins
    1173         792 :       IF (l < endbead) THEN
    1174             :          ! startbead belongs to the same atom
    1175         619 :          startatom = endatom
    1176         619 :          startbead = endbead - l
    1177             :       ELSE
    1178             :          ! startbead belongs to a different atom
    1179             :          ! find previous atom (assuming l < nbeads)
    1180         173 :          startatom = helium%iperm(endatom)
    1181         173 :          startbead = endbead + helium%beads - l
    1182             :       END IF
    1183             :       sold = worm_path_action(helium, helium%pos, &
    1184         792 :                               startatom, startbead, endatom, endbead)
    1185             : 
    1186         792 :       IF (helium%solute_present) THEN
    1187             :          ! yes this is correct, as the bead, that splits into tail and head only changes half
    1188             :          ! therefore only half of its action needs to be considred
    1189             :          ! and is cheated in here by passing it as head bead
    1190             :          sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
    1191             :                                                    startatom, startbead, &
    1192         767 :                                                    helium%pos(:, endatom, endbead), endatom, endbead)
    1193             :       END IF
    1194             : 
    1195         792 :       helium%worm_is_closed = .FALSE.
    1196         792 :       helium%worm_atom_idx_work = endatom
    1197         792 :       helium%worm_bead_idx_work = endbead
    1198             : 
    1199             :       ! alternative grow with consecutive gaussians
    1200         792 :       IF (startbead < endbead) THEN
    1201             :          ! everything belongs to the same atom
    1202             :          ! gro head from startbead
    1203        2123 :          DO kbead = startbead + 1, endbead - 1
    1204        6635 :             DO idim = 1, 3
    1205        4512 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1206        6016 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
    1207             :             END DO
    1208             :          END DO
    1209             :          ! last grow head bead
    1210        2476 :          DO idim = 1, 3
    1211        1857 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1212        2476 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
    1213             :          END DO
    1214         173 :       ELSE IF (endbead /= 1) THEN
    1215             :          ! is distributed among two atoms
    1216             :          ! grow from startbead
    1217         269 :          DO kbead = startbead + 1, helium%beads
    1218         692 :             DO idim = 1, 3
    1219         423 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1220         564 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
    1221             :             END DO
    1222             :          END DO
    1223             :          ! bead one of endatom relative to last on startatom
    1224         512 :          DO idim = 1, 3
    1225         384 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1226         512 :             helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
    1227             :          END DO
    1228             :          ! everything on endatom
    1229         243 :          DO kbead = 2, endbead - 1
    1230         588 :             DO idim = 1, 3
    1231         345 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1232         460 :                helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
    1233             :             END DO
    1234             :          END DO
    1235         512 :          DO idim = 1, 3
    1236         384 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1237         512 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, endatom, endbead - 1) + xr
    1238             :          END DO
    1239             :       ELSE ! imagtimewrap and headbead = 1
    1240             :          ! is distributed among two atoms
    1241             :          ! grow from startbead
    1242         165 :          DO kbead = startbead + 1, helium%beads
    1243         525 :             DO idim = 1, 3
    1244         360 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1245         480 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
    1246             :             END DO
    1247             :          END DO
    1248             :          ! bead one of endatom relative to last on startatom
    1249         180 :          DO idim = 1, 3
    1250         135 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1251         180 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, helium%beads) + xr
    1252             :          END DO
    1253             :       END IF
    1254             : 
    1255             :       snew = worm_path_action_worm_corrected(helium, helium%work, &
    1256             :                                              startatom, startbead, endatom, endbead, &
    1257         792 :                                              helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    1258             : 
    1259         792 :       IF (helium%solute_present) THEN
    1260             :          snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
    1261             :                                                    startatom, startbead, &
    1262         767 :                                                    helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    1263             :       END IF
    1264             : 
    1265             :       ! Metropolis:
    1266             :       ! first compute ln of free density matrix
    1267        3168 :       distvec(:) = helium%pos(:, startatom, startbead) - helium%pos(:, endatom, endbead)
    1268         792 :       CALL helium_pbc(helium, distvec)
    1269        3168 :       distsq = DOT_PRODUCT(distvec, distvec)
    1270             :       ! action difference
    1271         792 :       sdiff = sold - snew
    1272             :       ! modify action difference due to extra bead
    1273         792 :       sdiff = sdiff + distsq/(2.0_dp*helium%hb2m*REAL(l, dp)*helium%tau)
    1274         792 :       sdiff = sdiff + 1.5_dp*LOG(REAL(l, dp)*helium%tau)
    1275         792 :       sdiff = sdiff + helium%worm_ln_openclose_scale
    1276         792 :       IF (sdiff < 0) THEN
    1277         555 :          should_reject = .FALSE.
    1278         555 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    1279             :             should_reject = .TRUE.
    1280             :          ELSE
    1281         555 :             rtmp = helium%rng_stream_uniform%next()
    1282         555 :             IF (EXP(sdiff) < rtmp) THEN
    1283             :                should_reject = .TRUE.
    1284             :             END IF
    1285             :          END IF
    1286             :          IF (should_reject) THEN
    1287             :             ! rejected !
    1288             :             ! write back only changed atoms
    1289             :             ! transfer the new coordinates to work array
    1290         280 :             IF (startbead < endbead) THEN
    1291             :                ! everything belongs to the same atom
    1292         657 :                DO kbead = startbead + 1, endbead - 1
    1293        1953 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1294             :                END DO
    1295             :             ELSE
    1296             :                ! is distributed among two atoms
    1297             :                ! transfer to atom not containing the head bead
    1298         125 :                DO kbead = startbead + 1, helium%beads
    1299         335 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    1300             :                END DO
    1301             :                ! transfer to atom containing the head bead
    1302         123 :                DO kbead = 1, endbead - 1
    1303         327 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1304             :                END DO
    1305             :             END IF
    1306         280 :             helium%worm_is_closed = .TRUE.
    1307         280 :             ac = 0
    1308         283 :             RETURN
    1309             :          END IF
    1310             :       END IF
    1311             : 
    1312             :       ! for now accepted
    1313             :       ! rejection of the whole move if any centroid is farther away
    1314             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    1315             :       ! AND ist not moving towards the center
    1316         512 :       IF (.NOT. helium%periodic) THEN
    1317         498 :          IF (helium%solute_present) THEN
    1318        1992 :             new_com(:) = helium%center(:)
    1319         498 :             old_com(:) = new_com(:)
    1320             :          ELSE
    1321           0 :             new_com(:) = 0.0_dp
    1322           0 :             DO ia = 1, helium%atoms
    1323           0 :                DO ib = 1, helium%beads
    1324           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    1325             :                END DO
    1326             :             END DO
    1327           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1328             :             ! also compute the old center of mass (optimization potential)
    1329           0 :             old_com(:) = 0.0_dp
    1330           0 :             DO ia = 1, helium%atoms
    1331           0 :                DO ib = 1, helium%beads
    1332           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    1333             :                END DO
    1334             :             END DO
    1335           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1336             :          END IF
    1337         498 :          should_reject = .FALSE.
    1338       14772 :          atom: DO ia = 1, helium%atoms
    1339       14277 :             dr(:) = 0.0_dp
    1340      242709 :             DO ib = 1, helium%beads
    1341      928005 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    1342             :             END DO
    1343       57108 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    1344       57108 :             rtmp = DOT_PRODUCT(dr, dr)
    1345       14772 :             IF (rtmp >= helium%droplet_radius**2) THEN
    1346           3 :                dro(:) = 0.0_dp
    1347          51 :                DO ib = 1, helium%beads
    1348         195 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    1349             :                END DO
    1350          12 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    1351          12 :                rtmpo = DOT_PRODUCT(dro, dro)
    1352             :                ! only reject if it moves away from COM outside the droplet radius
    1353           3 :                IF (rtmpo < rtmp) THEN
    1354             :                   ! found - this one does not fit within R from the center
    1355             :                   should_reject = .TRUE.
    1356             :                   EXIT atom
    1357             :                END IF
    1358             :             END IF
    1359             :          END DO atom
    1360         498 :          IF (should_reject) THEN
    1361             :             ! restore original coordinates
    1362             :             ! write back only changed atoms
    1363           3 :             IF (startbead < endbead) THEN
    1364             :                ! everything belongs to the same atom
    1365           0 :                DO kbead = startbead + 1, endbead - 1
    1366           0 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1367             :                END DO
    1368             :             ELSE
    1369             :                ! is distributed among two atoms
    1370             :                ! transfer to atom not containing the head bead
    1371           6 :                DO kbead = startbead + 1, helium%beads
    1372          15 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    1373             :                END DO
    1374             :                ! transfer to atom containing the head bead
    1375           9 :                DO kbead = 1, endbead - 1
    1376          27 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1377             :                END DO
    1378             :             END IF
    1379           3 :             helium%worm_is_closed = .TRUE.
    1380             :             !helium%worm_atom_idx_work = helium%worm_atom_idx
    1381             :             !helium%worm_bead_idx_work = helium%worm_bead_idx
    1382           3 :             ac = 0
    1383           3 :             RETURN
    1384             :          END IF
    1385             :       END IF
    1386             : 
    1387             :       ! finally accepted
    1388         509 :       ac = 1
    1389             :       ! write changed coordinates to position array
    1390         509 :       IF (startbead < endbead) THEN
    1391             :          ! everything belongs to the same atom
    1392        1466 :          DO kbead = startbead + 1, endbead - 1
    1393        4682 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    1394             :          END DO
    1395             :       ELSE
    1396             :          ! is distributed among two atoms
    1397             :          ! transfer to atom not containing the head bead
    1398         303 :          DO kbead = startbead + 1, helium%beads
    1399         867 :             helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
    1400             :          END DO
    1401             :          ! transfer to atom containing the head bead
    1402         284 :          DO kbead = 1, endbead - 1
    1403         791 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    1404             :          END DO
    1405             :       END IF
    1406        2036 :       helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
    1407         509 :       helium%worm_atom_idx = helium%worm_atom_idx_work
    1408         509 :       helium%worm_bead_idx = helium%worm_bead_idx_work
    1409             : 
    1410             :    END SUBROUTINE worm_open_move
    1411             : 
    1412             : ! **************************************************************************************************
    1413             : !> \brief Close move back to diagonal elements of density matrix, permutation fixed in closed state
    1414             : !> \param pint_env ...
    1415             : !> \param helium ...
    1416             : !> \param l ...
    1417             : !> \param ac ...
    1418             : !> \author fuhl
    1419             : ! **************************************************************************************************
    1420        2069 :    SUBROUTINE worm_close_move(pint_env, helium, l, ac)
    1421             : 
    1422             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    1423             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1424             :       INTEGER, INTENT(IN)                                :: l
    1425             :       INTEGER, INTENT(OUT)                               :: ac
    1426             : 
    1427             :       INTEGER                                            :: endatom, endbead, ia, ib, jbead, kbead, &
    1428             :                                                             startatom, startbead
    1429             :       LOGICAL                                            :: should_reject
    1430             :       REAL(KIND=dp)                                      :: distsq, mass, rtmp, rtmpo, sdiff, snew, &
    1431             :                                                             sold
    1432             :       REAL(KIND=dp), DIMENSION(3)                        :: distvec, dr, dro, new_com, old_com
    1433        2069 :       REAL(KIND=dp), DIMENSION(3, l-1)                   :: newsection
    1434             : 
    1435        2069 :       mass = helium%he_mass_au
    1436             : 
    1437        2069 :       endatom = helium%worm_atom_idx
    1438        2069 :       endbead = helium%worm_bead_idx
    1439             :       ! get index of the atom and bead, where the resampling of the head begins
    1440        2069 :       IF (l < endbead) THEN
    1441             :          ! startbead belongs to the same atom
    1442        1545 :          startatom = endatom
    1443        1545 :          startbead = endbead - l
    1444             :       ELSE
    1445             :          ! startbead belongs to a different atom
    1446             :          ! find previous atom (assuming l < nbeads)
    1447         524 :          startatom = helium%iperm(endatom)
    1448         524 :          startbead = endbead + helium%beads - l
    1449             :       END IF
    1450             :       sold = worm_path_action_worm_corrected(helium, helium%pos, &
    1451             :                                              startatom, startbead, endatom, endbead, &
    1452        2069 :                                              helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
    1453             : 
    1454        2069 :       IF (helium%solute_present) THEN
    1455             :          sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
    1456             :                                                    startatom, startbead, &
    1457        2039 :                                                    helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
    1458             :       END IF
    1459             : 
    1460             :       ! close between head and tail
    1461             :       ! only l-1 beads need to be reconstructed
    1462             :       CALL path_construct(helium, &
    1463             :                           helium%pos(:, startatom, startbead), &
    1464             :                           helium%pos(:, endatom, endbead), l - 1, &
    1465        2069 :                           newsection)
    1466             : 
    1467             :       ! transfer the new coordinates to work array
    1468        2069 :       jbead = 1
    1469        2069 :       IF (startbead < endbead) THEN
    1470             :          ! everything belongs to the same atom
    1471        5304 :          DO kbead = startbead + 1, endbead - 1
    1472       15036 :             helium%work(:, endatom, kbead) = newsection(:, jbead)
    1473        5304 :             jbead = jbead + 1
    1474             :          END DO
    1475             :       ELSE
    1476             :          ! is distributed among two atoms
    1477             :          ! transfer to atom not containing the head bead
    1478        1225 :          DO kbead = startbead + 1, helium%beads
    1479        2804 :             helium%work(:, startatom, kbead) = newsection(:, jbead)
    1480        1225 :             jbead = jbead + 1
    1481             :          END DO
    1482             :          ! transfer to atom containing the head bead
    1483        1276 :          DO kbead = 1, endbead - 1
    1484        3008 :             helium%work(:, endatom, kbead) = newsection(:, jbead)
    1485        1276 :             jbead = jbead + 1
    1486             :          END DO
    1487             :       END IF
    1488             : 
    1489        2069 :       helium%worm_is_closed = .TRUE.
    1490             : 
    1491             :       snew = worm_path_action(helium, helium%work, &
    1492        2069 :                               startatom, startbead, endatom, endbead)
    1493             : 
    1494        2069 :       IF (helium%solute_present) THEN
    1495             :          ! yes this is correct, as the bead, that was split into tail and head only changes half
    1496             :          ! therefore only half of its action needs to be considred
    1497             :          ! and is cheated in here by passing it as head bead
    1498             :          snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
    1499             :                                                    startatom, startbead, &
    1500        2039 :                                                    helium%work(:, endatom, endbead), endatom, endbead)
    1501             :       END IF
    1502             : 
    1503             :       ! Metropolis:
    1504             :       ! first compute ln of free density matrix
    1505        8276 :       distvec(:) = helium%pos(:, startatom, startbead) - helium%pos(:, endatom, endbead)
    1506        2069 :       CALL helium_pbc(helium, distvec)
    1507        8276 :       distsq = DOT_PRODUCT(distvec, distvec)
    1508             :       ! action difference
    1509        2069 :       sdiff = sold - snew
    1510             :       ! modify action difference due to extra bead
    1511        2069 :       sdiff = sdiff - distsq/(2.0_dp*helium%hb2m*REAL(l, dp)*helium%tau)
    1512        2069 :       sdiff = sdiff - 1.5_dp*LOG(REAL(l, dp)*helium%tau)
    1513        2069 :       sdiff = sdiff - helium%worm_ln_openclose_scale
    1514        2069 :       IF (sdiff < 0) THEN
    1515        1834 :          should_reject = .FALSE.
    1516        1834 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    1517             :             should_reject = .TRUE.
    1518             :          ELSE
    1519        1798 :             rtmp = helium%rng_stream_uniform%next()
    1520        1798 :             IF (EXP(sdiff) < rtmp) THEN
    1521             :                should_reject = .TRUE.
    1522             :             END IF
    1523             :          END IF
    1524             :          IF (should_reject) THEN
    1525             :             ! rejected !
    1526             :             ! write back only changed atoms
    1527             :             ! transfer the new coordinates to work array
    1528        1559 :             IF (startbead < endbead) THEN
    1529             :                ! everything belongs to the same atom
    1530        3881 :                DO kbead = startbead + 1, endbead - 1
    1531       12056 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1532             :                END DO
    1533             :             ELSE
    1534             :                ! is distributed among two atoms
    1535             :                ! transfer to atom not containing the head bead
    1536         916 :                DO kbead = startbead + 1, helium%beads
    1537        2455 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    1538             :                END DO
    1539             :                ! transfer to atom containing the head bead
    1540         985 :                DO kbead = 1, endbead - 1
    1541        2731 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1542             :                END DO
    1543             :             END IF
    1544        1559 :             helium%worm_is_closed = .FALSE.
    1545        1559 :             ac = 0
    1546        1559 :             RETURN
    1547             :          END IF
    1548             :       END IF
    1549             : 
    1550             :       ! for now accepted
    1551             :       ! rejection of the whole move if any centroid is farther away
    1552             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    1553             :       ! AND ist not moving towards the center
    1554         510 :       IF (.NOT. helium%periodic) THEN
    1555         496 :          IF (helium%solute_present) THEN
    1556        1984 :             new_com(:) = helium%center(:)
    1557         496 :             old_com(:) = new_com(:)
    1558             :          ELSE
    1559           0 :             new_com(:) = 0.0_dp
    1560           0 :             DO ia = 1, helium%atoms
    1561           0 :                DO ib = 1, helium%beads
    1562           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    1563             :                END DO
    1564             :             END DO
    1565           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1566             :             ! also compute the old center of mass (optimization potential)
    1567           0 :             old_com(:) = 0.0_dp
    1568           0 :             DO ia = 1, helium%atoms
    1569           0 :                DO ib = 1, helium%beads
    1570           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    1571             :                END DO
    1572             :             END DO
    1573           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1574             :          END IF
    1575         496 :          should_reject = .FALSE.
    1576       14746 :          atom: DO ia = 1, helium%atoms
    1577       14251 :             dr(:) = 0.0_dp
    1578      242267 :             DO ib = 1, helium%beads
    1579      926315 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    1580             :             END DO
    1581       57004 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    1582       57004 :             rtmp = DOT_PRODUCT(dr, dr)
    1583       14746 :             IF (rtmp >= helium%droplet_radius**2) THEN
    1584           1 :                dro(:) = 0.0_dp
    1585          17 :                DO ib = 1, helium%beads
    1586          65 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    1587             :                END DO
    1588           4 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    1589           4 :                rtmpo = DOT_PRODUCT(dro, dro)
    1590             :                ! only reject if it moves away from COM outside the droplet radius
    1591           1 :                IF (rtmpo < rtmp) THEN
    1592             :                   ! found - this one does not fit within R from the center
    1593             :                   should_reject = .TRUE.
    1594             :                   EXIT atom
    1595             :                END IF
    1596             :             END IF
    1597             :          END DO atom
    1598         496 :          IF (should_reject) THEN
    1599             :             ! restore original coordinates
    1600             :             ! write back only changed atoms
    1601           1 :             IF (startbead < endbead) THEN
    1602             :                ! everything belongs to the same atom
    1603           5 :                DO kbead = startbead + 1, endbead - 1
    1604          17 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1605             :                END DO
    1606             :             ELSE
    1607             :                ! is distributed among two atoms
    1608             :                ! transfer to atom not containing the head bead
    1609           0 :                DO kbead = startbead + 1, helium%beads
    1610           0 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    1611             :                END DO
    1612             :                ! transfer to atom containing the head bead
    1613           0 :                DO kbead = 1, endbead - 1
    1614           0 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1615             :                END DO
    1616             :             END IF
    1617           1 :             helium%worm_is_closed = .FALSE.
    1618           1 :             ac = 0
    1619           1 :             RETURN
    1620             :          END IF
    1621             :       END IF
    1622             : 
    1623             :       ! finally accepted
    1624         509 :       ac = 1
    1625             :       ! write changed coordinates to position array
    1626         509 :       IF (startbead < endbead) THEN
    1627             :          ! everything belongs to the same atom
    1628        1418 :          DO kbead = startbead + 1, endbead - 1
    1629        4508 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    1630             :          END DO
    1631             :       ELSE
    1632             :          ! is distributed among two atoms
    1633             :          ! transfer to atom not containing the head bead
    1634         309 :          DO kbead = startbead + 1, helium%beads
    1635         873 :             helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
    1636             :          END DO
    1637             :          ! transfer to atom containing the head bead
    1638         291 :          DO kbead = 1, endbead - 1
    1639         801 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    1640             :          END DO
    1641             :       END IF
    1642             : 
    1643             :    END SUBROUTINE worm_close_move
    1644             : 
    1645             : ! **************************************************************************************************
    1646             : !> \brief Move head in open state
    1647             : !> \param pint_env ...
    1648             : !> \param helium ...
    1649             : !> \param l ...
    1650             : !> \param ac ...
    1651             : !> \author fuhl
    1652             : ! **************************************************************************************************
    1653        1033 :    SUBROUTINE worm_head_move(pint_env, helium, l, ac)
    1654             : 
    1655             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    1656             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1657             :       INTEGER, INTENT(IN)                                :: l
    1658             :       INTEGER, INTENT(OUT)                               :: ac
    1659             : 
    1660             :       INTEGER                                            :: endatom, endbead, ia, ib, idim, kbead, &
    1661             :                                                             startatom, startbead
    1662             :       LOGICAL                                            :: should_reject
    1663             :       REAL(KIND=dp)                                      :: mass, rtmp, rtmpo, sdiff, snew, sold, xr
    1664             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
    1665             : 
    1666        1033 :       mass = helium%he_mass_au
    1667             : 
    1668             :       ! get index of the atom and bead, where the resampling of the head begins
    1669        1033 :       endatom = helium%worm_atom_idx
    1670        1033 :       endbead = helium%worm_bead_idx
    1671        1033 :       IF (l < endbead) THEN
    1672             :          ! startbead belongs to the same atom
    1673         806 :          startatom = endatom
    1674         806 :          startbead = endbead - l
    1675             :       ELSE
    1676             :          ! startbead belongs to a different atom
    1677             :          ! find previous atom (assuming l < nbeads)
    1678         227 :          startatom = helium%iperm(endatom)
    1679         227 :          startbead = endbead + helium%beads - l
    1680             :       END IF
    1681             : 
    1682             :       sold = worm_path_action_worm_corrected(helium, helium%pos, &
    1683             :                                              startatom, startbead, endatom, endbead, &
    1684        1033 :                                              helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
    1685             : 
    1686        1033 :       IF (helium%solute_present) THEN
    1687             :          sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
    1688             :                                                    startatom, startbead, &
    1689        1021 :                                                    helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
    1690             :       END IF
    1691             : 
    1692             :       ! alternative grow with consecutive gaussians
    1693        1033 :       IF (startbead < endbead) THEN
    1694             :          ! everything belongs to the same atom
    1695             :          ! gro head from startbead
    1696        2751 :          DO kbead = startbead + 1, endbead - 1
    1697        8586 :             DO idim = 1, 3
    1698        5835 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1699        7780 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
    1700             :             END DO
    1701             :          END DO
    1702             :          ! last grow head bead
    1703        3224 :          DO idim = 1, 3
    1704        2418 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1705        3224 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
    1706             :          END DO
    1707         227 :       ELSE IF (endbead /= 1) THEN
    1708             :          ! is distributed among two atoms
    1709             :          ! grow from startbead
    1710         310 :          DO kbead = startbead + 1, helium%beads
    1711         781 :             DO idim = 1, 3
    1712         471 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1713         628 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
    1714             :             END DO
    1715             :          END DO
    1716             :          ! bead one of endatom relative to last on startatom
    1717         612 :          DO idim = 1, 3
    1718         459 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1719         612 :             helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
    1720             :          END DO
    1721             :          ! everything on endatom
    1722         300 :          DO kbead = 2, endbead - 1
    1723         741 :             DO idim = 1, 3
    1724         441 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1725         588 :                helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
    1726             :             END DO
    1727             :          END DO
    1728         612 :          DO idim = 1, 3
    1729         459 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1730         612 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, endatom, endbead - 1) + xr
    1731             :          END DO
    1732             :       ELSE ! imagtimewrap and headbead = 1
    1733             :          ! is distributed among two atoms
    1734             :          ! grow from startbead
    1735         252 :          DO kbead = startbead + 1, helium%beads
    1736         786 :             DO idim = 1, 3
    1737         534 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1738         712 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
    1739             :             END DO
    1740             :          END DO
    1741             :          ! bead one of endatom relative to last on startatom
    1742         296 :          DO idim = 1, 3
    1743         222 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1744         296 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, helium%beads) + xr
    1745             :          END DO
    1746             :       END IF
    1747             : 
    1748             :       snew = worm_path_action_worm_corrected(helium, helium%work, &
    1749             :                                              startatom, startbead, endatom, endbead, &
    1750        1033 :                                              helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    1751             : 
    1752        1033 :       IF (helium%solute_present) THEN
    1753             :          snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
    1754             :                                                    startatom, startbead, &
    1755        1021 :                                                    helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    1756             :       END IF
    1757             : 
    1758             :       ! Metropolis:
    1759             :       ! action difference
    1760        1033 :       sdiff = sold - snew
    1761             :       ! modify action difference due to extra bead
    1762        1033 :       IF (sdiff < 0) THEN
    1763         568 :          should_reject = .FALSE.
    1764         568 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    1765             :             should_reject = .TRUE.
    1766             :          ELSE
    1767         561 :             rtmp = helium%rng_stream_uniform%next()
    1768         561 :             IF (EXP(sdiff) < rtmp) THEN
    1769             :                should_reject = .TRUE.
    1770             :             END IF
    1771             :          END IF
    1772             :          IF (should_reject) THEN
    1773             :             ! rejected !
    1774             :             ! write back only changed atoms
    1775             :             ! transfer the new coordinates to work array
    1776          32 :             IF (startbead < endbead) THEN
    1777             :                ! everything belongs to the same atom
    1778         110 :                DO kbead = startbead + 1, endbead - 1
    1779         353 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1780             :                END DO
    1781             :             ELSE
    1782             :                ! is distributed among two atoms
    1783             :                ! transfer to atom not containing the head bead
    1784           9 :                DO kbead = startbead + 1, helium%beads
    1785          27 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    1786             :                END DO
    1787             :                ! transfer to atom containing the head bead
    1788           9 :                DO kbead = 1, endbead - 1
    1789          27 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1790             :                END DO
    1791             :             END IF
    1792         128 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    1793          32 :             ac = 0
    1794          36 :             RETURN
    1795             :          END IF
    1796             :       END IF
    1797             : 
    1798             :       ! for now accepted
    1799             :       ! rejection of the whole move if any centroid is farther away
    1800             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    1801             :       ! AND ist not moving towards the center
    1802        1001 :       IF (.NOT. helium%periodic) THEN
    1803         994 :          IF (helium%solute_present) THEN
    1804        3976 :             new_com(:) = helium%center(:)
    1805         994 :             old_com(:) = new_com(:)
    1806             :          ELSE
    1807           0 :             new_com(:) = 0.0_dp
    1808           0 :             DO ia = 1, helium%atoms
    1809           0 :                DO ib = 1, helium%beads
    1810           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    1811             :                END DO
    1812             :             END DO
    1813           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1814             :             ! also compute the old center of mass (optimization potential)
    1815           0 :             old_com(:) = 0.0_dp
    1816           0 :             DO ia = 1, helium%atoms
    1817           0 :                DO ib = 1, helium%beads
    1818           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    1819             :                END DO
    1820             :             END DO
    1821           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    1822             :          END IF
    1823         994 :          should_reject = .FALSE.
    1824       30520 :          atom: DO ia = 1, helium%atoms
    1825       29530 :             dr(:) = 0.0_dp
    1826      502010 :             DO ib = 1, helium%beads
    1827     1919450 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    1828             :             END DO
    1829      118120 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    1830      118120 :             rtmp = DOT_PRODUCT(dr, dr)
    1831       30520 :             IF (rtmp >= helium%droplet_radius**2) THEN
    1832           4 :                dro(:) = 0.0_dp
    1833          68 :                DO ib = 1, helium%beads
    1834         260 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    1835             :                END DO
    1836          16 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    1837          16 :                rtmpo = DOT_PRODUCT(dro, dro)
    1838             :                ! only reject if it moves away from COM outside the droplet radius
    1839           4 :                IF (rtmpo < rtmp) THEN
    1840             :                   ! found - this one does not fit within R from the center
    1841             :                   should_reject = .TRUE.
    1842             :                   EXIT atom
    1843             :                END IF
    1844             :             END IF
    1845             :          END DO atom
    1846         994 :          IF (should_reject) THEN
    1847             :             ! restore original coordinates
    1848             :             ! write back only changed atoms
    1849           4 :             IF (startbead < endbead) THEN
    1850             :                ! everything belongs to the same atom
    1851          13 :                DO kbead = startbead + 1, endbead - 1
    1852          43 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1853             :                END DO
    1854             :             ELSE
    1855             :                ! is distributed among two atoms
    1856             :                ! transfer to atom not containing the head bead
    1857           1 :                DO kbead = startbead + 1, helium%beads
    1858           1 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    1859             :                END DO
    1860             :                ! transfer to atom containing the head bead
    1861           5 :                DO kbead = 1, endbead - 1
    1862          17 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    1863             :                END DO
    1864             :             END IF
    1865          16 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    1866           4 :             ac = 0
    1867           4 :             RETURN
    1868             :          END IF
    1869             :       END IF
    1870             : 
    1871             :       ! finally accepted
    1872         997 :       ac = 1
    1873             :       ! write changed coordinates to position array
    1874         997 :       IF (startbead < endbead) THEN
    1875             :          ! everything belongs to the same atom
    1876        2628 :          DO kbead = startbead + 1, endbead - 1
    1877        8190 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    1878             :          END DO
    1879             :       ELSE
    1880             :          ! is distributed among two atoms
    1881             :          ! transfer to atom not containing the head bead
    1882         552 :          DO kbead = startbead + 1, helium%beads
    1883        1539 :             helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
    1884             :          END DO
    1885             :          ! transfer to atom containing the head bead
    1886         513 :          DO kbead = 1, endbead - 1
    1887        1383 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    1888             :          END DO
    1889             :       END IF
    1890        3988 :       helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
    1891             : 
    1892             :    END SUBROUTINE worm_head_move
    1893             : 
    1894             : ! **************************************************************************************************
    1895             : !> \brief Move tail in open state
    1896             : !> \param pint_env ...
    1897             : !> \param helium ...
    1898             : !> \param l ...
    1899             : !> \param ac ...
    1900             : !> \author fuhl
    1901             : ! **************************************************************************************************
    1902        1041 :    SUBROUTINE worm_tail_move(pint_env, helium, l, ac)
    1903             : 
    1904             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    1905             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1906             :       INTEGER, INTENT(IN)                                :: l
    1907             :       INTEGER, INTENT(OUT)                               :: ac
    1908             : 
    1909             :       INTEGER                                            :: endatom, endbead, ia, ib, idim, kbead, &
    1910             :                                                             startatom, startbead
    1911             :       LOGICAL                                            :: should_reject
    1912             :       REAL(KIND=dp)                                      :: mass, rtmp, rtmpo, sdiff, snew, sold, xr
    1913             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
    1914             : 
    1915        1041 :       mass = helium%he_mass_au
    1916             : 
    1917             :       ! get index of the atom and bead, where the resampling of the tail ends
    1918        1041 :       startatom = helium%worm_atom_idx
    1919        1041 :       startbead = helium%worm_bead_idx
    1920        1041 :       endbead = startbead + l
    1921             : 
    1922        1041 :       IF (endbead <= helium%beads) THEN
    1923             :          ! endbead belongs to the same atom
    1924         862 :          endatom = startatom
    1925             :       ELSE
    1926             :          ! endbead belongs to a different atom
    1927             :          ! find next atom (assuming l < nbeads)
    1928         179 :          endatom = helium%permutation(startatom)
    1929         179 :          endbead = endbead - helium%beads
    1930             :       END IF
    1931             : 
    1932             :       !yes this is correct, as the head does not play any role here
    1933             :       sold = worm_path_action(helium, helium%pos, &
    1934        1041 :                               startatom, startbead, endatom, endbead)
    1935             : 
    1936        1041 :       IF (helium%solute_present) THEN
    1937             :          sold = sold + worm_path_inter_action_tail(pint_env, helium, helium%pos, &
    1938             :                                                    endatom, endbead, &
    1939        1029 :                                                    helium%worm_atom_idx, helium%worm_bead_idx)
    1940             :       END IF
    1941             : 
    1942             :       ! alternative grow with consecutive gaussians
    1943        1041 :       IF (startbead < endbead) THEN
    1944             :          ! everything belongs to the same atom
    1945             :          ! gro tail from endbead to startbead (confusing eh?)
    1946        3777 :          DO kbead = endbead - 1, startbead, -1
    1947       12522 :             DO idim = 1, 3
    1948        8745 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1949       11660 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead + 1) + xr
    1950             :             END DO
    1951             :          END DO
    1952             :       ELSE
    1953             :          ! is distributed among two atoms
    1954             :          ! grow from endbead
    1955         379 :          DO kbead = endbead - 1, 1, -1
    1956         979 :             DO idim = 1, 3
    1957         600 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1958         800 :                helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead + 1) + xr
    1959             :             END DO
    1960             :          END DO
    1961             : 
    1962             :          ! over imaginary time boundary
    1963         716 :          DO idim = 1, 3
    1964         537 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1965         716 :             helium%work(idim, startatom, helium%beads) = helium%work(idim, endatom, 1) + xr
    1966             :          END DO
    1967             : 
    1968             :          ! rest on startatom
    1969         480 :          DO kbead = helium%beads - 1, startbead, -1
    1970        1383 :             DO idim = 1, 3
    1971         903 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    1972        1204 :                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead + 1) + xr
    1973             :             END DO
    1974             :          END DO
    1975             :       END IF
    1976             : 
    1977             :       !yes this is correct, as the head does not play any role here
    1978             :       snew = worm_path_action(helium, helium%work, &
    1979        1041 :                               startatom, startbead, endatom, endbead)
    1980             : 
    1981        1041 :       IF (helium%solute_present) THEN
    1982             :          snew = snew + worm_path_inter_action_tail(pint_env, helium, helium%work, &
    1983             :                                                    endatom, endbead, &
    1984        1029 :                                                    helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    1985             :       END IF
    1986             : 
    1987             :       ! Metropolis:
    1988             :       ! action difference
    1989        1041 :       sdiff = sold - snew
    1990             :       ! modify action difference due to extra bead
    1991        1041 :       IF (sdiff < 0) THEN
    1992         526 :          should_reject = .FALSE.
    1993         526 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    1994             :             should_reject = .TRUE.
    1995             :          ELSE
    1996         526 :             rtmp = helium%rng_stream_uniform%next()
    1997         526 :             IF (EXP(sdiff) < rtmp) THEN
    1998             :                should_reject = .TRUE.
    1999             :             END IF
    2000             :          END IF
    2001             :          IF (should_reject) THEN
    2002             :             ! rejected !
    2003             :             ! write back only changed atoms
    2004             :             ! transfer the new coordinates to work array
    2005          25 :             IF (startbead < endbead) THEN
    2006             :                ! everything belongs to the same atom
    2007         100 :                DO kbead = startbead, endbead - 1
    2008         343 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    2009             :                END DO
    2010             :             ELSE
    2011             :                ! is distributed among two atoms
    2012             :                ! transfer to atom not containing the tail bead
    2013          32 :                DO kbead = startbead, helium%beads
    2014         110 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    2015             :                END DO
    2016             :                ! transfer to atom containing the tail bead
    2017           8 :                DO kbead = 1, endbead - 1
    2018          14 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    2019             :                END DO
    2020             :             END IF
    2021         100 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2022          25 :             ac = 0
    2023          32 :             RETURN
    2024             :          END IF
    2025             :       END IF
    2026             : 
    2027             :       ! for now accepted
    2028             :       ! rejection of the whole move if any centroid is farther away
    2029             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    2030             :       ! AND ist not moving towards the center
    2031        1016 :       IF (.NOT. helium%periodic) THEN
    2032        1010 :          IF (helium%solute_present) THEN
    2033        4040 :             new_com(:) = helium%center(:)
    2034        1010 :             old_com(:) = new_com(:)
    2035             :          ELSE
    2036           0 :             new_com(:) = 0.0_dp
    2037           0 :             DO ia = 1, helium%atoms
    2038           0 :                DO ib = 1, helium%beads
    2039           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    2040             :                END DO
    2041             :             END DO
    2042           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2043             :             ! also compute the old center of mass (optimization potential)
    2044           0 :             old_com(:) = 0.0_dp
    2045           0 :             DO ia = 1, helium%atoms
    2046           0 :                DO ib = 1, helium%beads
    2047           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    2048             :                END DO
    2049             :             END DO
    2050           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2051             :          END IF
    2052        1010 :          should_reject = .FALSE.
    2053       31174 :          atom: DO ia = 1, helium%atoms
    2054       30171 :             dr(:) = 0.0_dp
    2055      512907 :             DO ib = 1, helium%beads
    2056     1961115 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    2057             :             END DO
    2058      120684 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    2059      120684 :             rtmp = DOT_PRODUCT(dr, dr)
    2060       31174 :             IF (rtmp >= helium%droplet_radius**2) THEN
    2061           7 :                dro(:) = 0.0_dp
    2062         119 :                DO ib = 1, helium%beads
    2063         455 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    2064             :                END DO
    2065          28 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    2066          28 :                rtmpo = DOT_PRODUCT(dro, dro)
    2067             :                ! only reject if it moves away from COM outside the droplet radius
    2068           7 :                IF (rtmpo < rtmp) THEN
    2069             :                   ! found - this one does not fit within R from the center
    2070             :                   should_reject = .TRUE.
    2071             :                   EXIT atom
    2072             :                END IF
    2073             :             END IF
    2074             :          END DO atom
    2075        1010 :          IF (should_reject) THEN
    2076             :             ! restore original coordinates
    2077             :             ! write back only changed atoms
    2078           7 :             IF (startbead < endbead) THEN
    2079             :                ! everything belongs to the same atom
    2080          30 :                DO kbead = startbead, endbead - 1
    2081          99 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    2082             :                END DO
    2083             :             ELSE
    2084             :                ! is distributed among two atoms
    2085             :                ! transfer to atom not containing the tail bead
    2086           0 :                DO kbead = startbead, helium%beads
    2087           0 :                   helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
    2088             :                END DO
    2089             :                ! transfer to atom containing the tail bead
    2090           0 :                DO kbead = 1, endbead - 1
    2091           0 :                   helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
    2092             :                END DO
    2093             :             END IF
    2094          28 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2095           7 :             ac = 0
    2096           7 :             RETURN
    2097             :          END IF
    2098             :       END IF
    2099             : 
    2100             :       ! finally accepted
    2101        1009 :       ac = 1
    2102             :       ! write changed coordinates to position array
    2103        1009 :       IF (startbead < endbead) THEN
    2104             :          ! everything belongs to the same atom
    2105        3647 :          DO kbead = startbead, endbead - 1
    2106       12080 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    2107             :          END DO
    2108             :       ELSE
    2109             :          ! is distributed among two atoms
    2110             :          ! transfer to atom not containing the tail bead
    2111         627 :          DO kbead = startbead, helium%beads
    2112        1989 :             helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
    2113             :          END DO
    2114             :          ! transfer to atom containing the tail bead
    2115         371 :          DO kbead = 1, endbead - 1
    2116         965 :             helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
    2117             :          END DO
    2118             :       END IF
    2119             : 
    2120             :    END SUBROUTINE worm_tail_move
    2121             : 
    2122             : ! **************************************************************************************************
    2123             : !> \brief Crawl forward in open state, can avoid being trapped in open state
    2124             : !> \param pint_env ...
    2125             : !> \param helium ...
    2126             : !> \param l ...
    2127             : !> \param ac ...
    2128             : !> \author fuhl
    2129             : ! **************************************************************************************************
    2130        1880 :    SUBROUTINE worm_crawl_move_forward(pint_env, helium, l, ac)
    2131             : 
    2132             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2133             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    2134             :       INTEGER, INTENT(IN)                                :: l
    2135             :       INTEGER, INTENT(OUT)                               :: ac
    2136             : 
    2137             :       INTEGER                                            :: ia, ib, idim, kbead
    2138             :       LOGICAL                                            :: should_reject
    2139             :       REAL(KIND=dp)                                      :: mass, newtailpot, oldheadpot, &
    2140             :                                                             oldtailpot, rtmp, rtmpo, sdiff, snew, &
    2141             :                                                             sold, xr
    2142             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
    2143             : 
    2144        1880 :       mass = helium%he_mass_au
    2145             : 
    2146             :       ! determine position of new head in imaginary time
    2147        1880 :       helium%worm_bead_idx_work = helium%worm_bead_idx + l
    2148        1880 :       IF (helium%worm_bead_idx_work > helium%beads) THEN
    2149         418 :          helium%worm_bead_idx_work = helium%worm_bead_idx_work - helium%beads
    2150         418 :          helium%worm_atom_idx_work = helium%permutation(helium%worm_atom_idx)
    2151             :       ELSE
    2152        1462 :          helium%worm_atom_idx_work = helium%worm_atom_idx
    2153             :       END IF
    2154             : 
    2155             :       ! compute action before move
    2156             :       ! head is not involved in pair action
    2157             :       sold = worm_path_action(helium, helium%pos, &
    2158             :                               helium%worm_atom_idx, helium%worm_bead_idx, &
    2159        1880 :                               helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    2160        1880 :       IF (helium%solute_present) THEN
    2161             :          !this will leave out the old and new tail bead
    2162             :          ! due to efficiency reasons they are treated separately
    2163             :          sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
    2164             :                                               helium%worm_atom_idx, helium%worm_bead_idx, &
    2165        1848 :                                               helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    2166             : 
    2167             :          ! compute old/new head/tail interactions
    2168             :          ! old tail
    2169             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2170             :                                      helium%worm_atom_idx, helium%worm_bead_idx, &
    2171             :                                      helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), &
    2172        1848 :                                      energy=oldtailpot)
    2173        1848 :          oldtailpot = oldtailpot*helium%tau
    2174             : 
    2175             :          ! new tail
    2176             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2177             :                                      helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
    2178             :                                      helium%pos(:, helium%worm_atom_idx_work, helium%worm_bead_idx_work), &
    2179        1848 :                                      energy=newtailpot)
    2180        1848 :          newtailpot = newtailpot*helium%tau
    2181             : 
    2182             :          ! old head
    2183             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2184             :                                      helium%worm_atom_idx, helium%worm_bead_idx, &
    2185             :                                      helium%worm_xtra_bead, &
    2186        1848 :                                      energy=oldheadpot)
    2187        1848 :          oldheadpot = oldheadpot*helium%tau
    2188             : 
    2189             :          ! new head is not known yet
    2190             : 
    2191        1848 :          sold = sold + 0.5_dp*(oldtailpot + oldheadpot) + newtailpot
    2192             :       END IF
    2193             : 
    2194             :       ! copy over old head position to working array and grow from there
    2195        7520 :       helium%work(:, helium%worm_atom_idx, helium%worm_bead_idx) = helium%worm_xtra_bead
    2196             : 
    2197             :       ! grow head part with consecutive gaussians
    2198        1880 :       IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
    2199             :          ! everything belongs to the same atom
    2200             :          ! gro head from startbead
    2201        4999 :          DO kbead = helium%worm_bead_idx + 1, helium%worm_bead_idx_work - 1
    2202       15610 :             DO idim = 1, 3
    2203       10611 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2204       14148 :                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
    2205             :             END DO
    2206             :          END DO
    2207             :          ! last grow head bead
    2208        5848 :          DO idim = 1, 3
    2209        4386 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2210        5848 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%worm_bead_idx_work - 1) + xr
    2211             :          END DO
    2212         418 :       ELSE IF (helium%worm_bead_idx_work /= 1) THEN
    2213             :          ! is distributed among two atoms
    2214             :          ! grow from startbead
    2215         593 :          DO kbead = helium%worm_bead_idx + 1, helium%beads
    2216        1484 :             DO idim = 1, 3
    2217         891 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2218        1188 :                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
    2219             :             END DO
    2220             :          END DO
    2221             :          ! bead one of endatom relative to last on helium%worm_atom_idx
    2222        1184 :          DO idim = 1, 3
    2223         888 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2224        1184 :             helium%work(idim, helium%worm_atom_idx_work, 1) = helium%work(idim, helium%worm_atom_idx, helium%beads) + xr
    2225             :          END DO
    2226             :          ! everything on endatom
    2227         573 :          DO kbead = 2, helium%worm_bead_idx_work - 1
    2228        1404 :             DO idim = 1, 3
    2229         831 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2230        1108 :                helium%work(idim, helium%worm_atom_idx_work, kbead) = helium%work(idim, helium%worm_atom_idx_work, kbead - 1) + xr
    2231             :             END DO
    2232             :          END DO
    2233        1184 :          DO idim = 1, 3
    2234         888 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2235        1184 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx_work, helium%worm_bead_idx_work - 1) + xr
    2236             :          END DO
    2237             :       ELSE ! imagtimewrap and headbead = 1
    2238             :          ! is distributed among two atoms
    2239             :          ! grow from startbead
    2240         402 :          DO kbead = helium%worm_bead_idx + 1, helium%beads
    2241        1242 :             DO idim = 1, 3
    2242         840 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2243        1120 :                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
    2244             :             END DO
    2245             :          END DO
    2246             :          ! bead one of endatom relative to last on helium%worm_atom_idx
    2247         488 :          DO idim = 1, 3
    2248         366 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2249         488 :             helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%beads) + xr
    2250             :          END DO
    2251             :       END IF
    2252             : 
    2253             :       snew = worm_path_action_worm_corrected(helium, helium%work, &
    2254             :                                              helium%worm_atom_idx, helium%worm_bead_idx, &
    2255             :                                              helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
    2256        1880 :                                              helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    2257             : 
    2258        1880 :       IF (helium%solute_present) THEN
    2259             :          snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
    2260             :                                                    helium%worm_atom_idx, helium%worm_bead_idx, &
    2261        1848 :                                                    helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    2262        1848 :          snew = snew + 0.5_dp*newtailpot + oldheadpot
    2263             :       END IF
    2264             : 
    2265             :       ! Metropolis:
    2266             :       ! action difference
    2267        1880 :       sdiff = sold - snew
    2268             :       ! modify action difference due to extra bead
    2269        1880 :       IF (sdiff < 0) THEN
    2270         969 :          should_reject = .FALSE.
    2271         969 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    2272             :             should_reject = .TRUE.
    2273             :          ELSE
    2274         959 :             rtmp = helium%rng_stream_uniform%next()
    2275         959 :             IF (EXP(sdiff) < rtmp) THEN
    2276             :                should_reject = .TRUE.
    2277             :             END IF
    2278             :          END IF
    2279             :          IF (should_reject) THEN
    2280             :             ! rejected !
    2281             :             ! write back only changed atoms
    2282             :             ! transfer the new coordinates to work array
    2283          95 :             IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
    2284             :                ! everything belongs to the same atom
    2285         328 :                DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
    2286        1099 :                   helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
    2287             :                END DO
    2288             :             ELSE
    2289             :                ! is distributed among two atoms
    2290             :                ! transfer to atom not containing the head bead
    2291          68 :                DO kbead = helium%worm_bead_idx, helium%beads
    2292         200 :                   helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
    2293             :                END DO
    2294             :                ! transfer to atom containing the head bead
    2295          69 :                DO kbead = 1, helium%worm_bead_idx_work - 1
    2296         204 :                   helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
    2297             :                END DO
    2298             :             END IF
    2299         380 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2300          95 :             helium%worm_bead_idx_work = helium%worm_bead_idx
    2301          95 :             helium%worm_atom_idx_work = helium%worm_atom_idx
    2302          95 :             ac = 0
    2303         108 :             RETURN
    2304             :          END IF
    2305             :       END IF
    2306             : 
    2307             :       ! for now accepted
    2308             :       ! rejection of the whole move if any centroid is farther away
    2309             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    2310             :       ! AND ist not moving towards the center
    2311        1785 :       IF (.NOT. helium%periodic) THEN
    2312        1773 :          IF (helium%solute_present) THEN
    2313        7092 :             new_com(:) = helium%center(:)
    2314        1773 :             old_com(:) = new_com(:)
    2315             :          ELSE
    2316           0 :             new_com(:) = 0.0_dp
    2317           0 :             DO ia = 1, helium%atoms
    2318           0 :                DO ib = 1, helium%beads
    2319           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    2320             :                END DO
    2321             :             END DO
    2322           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2323             :             ! also compute the old center of mass (optimization potential)
    2324           0 :             old_com(:) = 0.0_dp
    2325           0 :             DO ia = 1, helium%atoms
    2326           0 :                DO ib = 1, helium%beads
    2327           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    2328             :                END DO
    2329             :             END DO
    2330           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2331             :          END IF
    2332        1773 :          should_reject = .FALSE.
    2333       53384 :          atom: DO ia = 1, helium%atoms
    2334       51624 :             dr(:) = 0.0_dp
    2335      877608 :             DO ib = 1, helium%beads
    2336     3355560 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    2337             :             END DO
    2338      206496 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    2339      206496 :             rtmp = DOT_PRODUCT(dr, dr)
    2340       53384 :             IF (rtmp >= helium%droplet_radius**2) THEN
    2341          13 :                dro(:) = 0.0_dp
    2342         221 :                DO ib = 1, helium%beads
    2343         845 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    2344             :                END DO
    2345          52 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    2346          52 :                rtmpo = DOT_PRODUCT(dro, dro)
    2347             :                ! only reject if it moves away from COM outside the droplet radius
    2348          13 :                IF (rtmpo < rtmp) THEN
    2349             :                   ! found - this one does not fit within R from the center
    2350             :                   should_reject = .TRUE.
    2351             :                   EXIT atom
    2352             :                END IF
    2353             :             END IF
    2354             :          END DO atom
    2355        1773 :          IF (should_reject) THEN
    2356             :             ! restore original coordinates
    2357             :             ! write back only changed atoms
    2358          13 :             IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
    2359             :                ! everything belongs to the same atom
    2360          40 :                DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
    2361         133 :                   helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
    2362             :                END DO
    2363             :             ELSE
    2364             :                ! is distributed among two atoms
    2365             :                ! transfer to atom not containing the head bead
    2366          13 :                DO kbead = helium%worm_bead_idx, helium%beads
    2367          40 :                   helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
    2368             :                END DO
    2369             :                ! transfer to atom containing the head bead
    2370          12 :                DO kbead = 1, helium%worm_bead_idx_work - 1
    2371          36 :                   helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
    2372             :                END DO
    2373             :             END IF
    2374          52 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2375          13 :             helium%worm_bead_idx_work = helium%worm_bead_idx
    2376          13 :             helium%worm_atom_idx_work = helium%worm_atom_idx
    2377          13 :             ac = 0
    2378          13 :             RETURN
    2379             :          END IF
    2380             :       END IF
    2381             : 
    2382             :       ! finally accepted
    2383        1772 :       ac = 1
    2384             :       ! write changed coordinates to position array
    2385        1772 :       IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
    2386             :          ! everything belongs to the same atom
    2387        6093 :          DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
    2388       20226 :             helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
    2389             :          END DO
    2390             :       ELSE
    2391             :          ! is distributed among two atoms
    2392             :          ! transfer to atom not containing the head bead
    2393        1332 :          DO kbead = helium%worm_bead_idx, helium%beads
    2394        4158 :             helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
    2395             :          END DO
    2396             :          ! transfer to atom containing the head bead
    2397         910 :          DO kbead = 1, helium%worm_bead_idx_work - 1
    2398        2470 :             helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
    2399             :          END DO
    2400             :       END IF
    2401        7088 :       helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
    2402        1772 :       helium%worm_bead_idx = helium%worm_bead_idx_work
    2403        1772 :       helium%worm_atom_idx = helium%worm_atom_idx_work
    2404             : 
    2405             :    END SUBROUTINE worm_crawl_move_forward
    2406             : 
    2407             : ! **************************************************************************************************
    2408             : !> \brief Crawl backwards in open state, can avoid being trapped in open state
    2409             : !> \param pint_env ...
    2410             : !> \param helium ...
    2411             : !> \param l ...
    2412             : !> \param ac ...
    2413             : !> \author fuhl
    2414             : ! **************************************************************************************************
    2415        2066 :    SUBROUTINE worm_crawl_move_backward(pint_env, helium, l, ac)
    2416             : 
    2417             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2418             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    2419             :       INTEGER, INTENT(IN)                                :: l
    2420             :       INTEGER, INTENT(OUT)                               :: ac
    2421             : 
    2422             :       INTEGER                                            :: ia, ib, idim, kbead
    2423             :       LOGICAL                                            :: should_reject
    2424             :       REAL(KIND=dp)                                      :: mass, newheadpot, oldheadpot, &
    2425             :                                                             oldtailpot, rtmp, rtmpo, sdiff, snew, &
    2426             :                                                             sold, xr
    2427             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
    2428             : 
    2429        2066 :       mass = helium%he_mass_au
    2430             : 
    2431             :       ! determine position of new head in imaginary time
    2432        2066 :       helium%worm_bead_idx_work = helium%worm_bead_idx - l
    2433        2066 :       IF (helium%worm_bead_idx_work < 1) THEN
    2434         447 :          helium%worm_bead_idx_work = helium%worm_bead_idx_work + helium%beads
    2435         447 :          helium%worm_atom_idx_work = helium%iperm(helium%worm_atom_idx)
    2436             :       ELSE
    2437        1619 :          helium%worm_atom_idx_work = helium%worm_atom_idx
    2438             :       END IF
    2439             : 
    2440             :       ! compute action before move
    2441             :       ! head is not involved in pair action
    2442             :       sold = worm_path_action_worm_corrected(helium, helium%work, &
    2443             :                                              helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
    2444             :                                              helium%worm_atom_idx, helium%worm_bead_idx, &
    2445        2066 :                                              helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
    2446             : 
    2447        2066 :       IF (helium%solute_present) THEN
    2448             :          !this will leave out the old and new tail bead
    2449             :          ! due to efficiency reasons they are treated separately
    2450             :          sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
    2451             :                                               helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
    2452        2032 :                                               helium%worm_atom_idx, helium%worm_bead_idx)
    2453             : 
    2454             :          ! compute old/new head/tail interactions
    2455             :          ! old tail
    2456             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2457             :                                      helium%worm_atom_idx, helium%worm_bead_idx, &
    2458             :                                      helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), &
    2459        2032 :                                      energy=oldtailpot)
    2460        2032 :          oldtailpot = oldtailpot*helium%tau
    2461             : 
    2462             :          ! old head
    2463             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2464             :                                      helium%worm_atom_idx, helium%worm_bead_idx, &
    2465             :                                      helium%worm_xtra_bead, &
    2466        2032 :                                      energy=oldheadpot)
    2467        2032 :          oldheadpot = oldheadpot*helium%tau
    2468             : 
    2469             :          ! new head
    2470             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2471             :                                      helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
    2472             :                                      helium%pos(:, helium%worm_atom_idx_work, helium%worm_bead_idx_work), &
    2473        2032 :                                      energy=newheadpot)
    2474        2032 :          newheadpot = newheadpot*helium%tau
    2475             : 
    2476             :          ! new tail not known yet
    2477             : 
    2478        2032 :          sold = sold + 0.5_dp*(oldtailpot + oldheadpot) + newheadpot
    2479             :       END IF
    2480             : 
    2481             :       ! copy position to the head bead
    2482        8264 :       helium%worm_xtra_bead_work = helium%pos(:, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    2483             : 
    2484             :       ! alternative grow with consecutive gaussians
    2485        2066 :       IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
    2486             :          ! everything belongs to the same atom
    2487             :          ! gro tail from endbead to startbead (confusing eh?)
    2488        7197 :          DO kbead = helium%worm_bead_idx - 1, helium%worm_bead_idx_work, -1
    2489       23931 :             DO idim = 1, 3
    2490       16734 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2491       22312 :                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead + 1) + xr
    2492             :             END DO
    2493             : 
    2494             :          END DO
    2495             :       ELSE
    2496             :          ! is distributed among two atoms
    2497             :          ! grow from endbead
    2498        1092 :          DO kbead = helium%worm_bead_idx - 1, 1, -1
    2499        3027 :             DO idim = 1, 3
    2500        1935 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2501        2580 :                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead + 1) + xr
    2502             :             END DO
    2503             :          END DO
    2504             : 
    2505             :          ! over imaginary time boundary
    2506        1788 :          DO idim = 1, 3
    2507        1341 :             xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2508        1788 :             helium%work(idim, helium%worm_atom_idx_work, helium%beads) = helium%work(idim, helium%worm_atom_idx, 1) + xr
    2509             :          END DO
    2510             : 
    2511             :          ! rest on startatom
    2512        1075 :          DO kbead = helium%beads - 1, helium%worm_bead_idx_work, -1
    2513        2959 :             DO idim = 1, 3
    2514        1884 :                xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
    2515        2512 :                helium%work(idim, helium%worm_atom_idx_work, kbead) = helium%work(idim, helium%worm_atom_idx_work, kbead + 1) + xr
    2516             :             END DO
    2517             :          END DO
    2518             :       END IF
    2519             : 
    2520             :       snew = worm_path_action(helium, helium%work, &
    2521             :                               helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
    2522        2066 :                               helium%worm_atom_idx, helium%worm_bead_idx)
    2523             : 
    2524        2066 :       IF (helium%solute_present) THEN
    2525             :          snew = snew + worm_path_inter_action_tail(pint_env, helium, helium%work, &
    2526             :                                                    helium%worm_atom_idx, helium%worm_bead_idx, &
    2527        2032 :                                                    helium%worm_atom_idx_work, helium%worm_bead_idx_work)
    2528        2032 :          snew = snew + 0.5_dp*newheadpot + oldtailpot
    2529             :       END IF
    2530             : 
    2531             :       ! Metropolis:
    2532             :       ! action difference
    2533        2066 :       sdiff = sold - snew
    2534             :       ! modify action difference due to extra bead
    2535        2066 :       IF (sdiff < 0) THEN
    2536        1097 :          should_reject = .FALSE.
    2537        1097 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    2538             :             should_reject = .TRUE.
    2539             :          ELSE
    2540        1096 :             rtmp = helium%rng_stream_uniform%next()
    2541        1096 :             IF (EXP(sdiff) < rtmp) THEN
    2542             :                should_reject = .TRUE.
    2543             :             END IF
    2544             :          END IF
    2545             :          IF (should_reject) THEN
    2546             :             ! rejected !
    2547             :             ! write back only changed atoms
    2548             :             ! transfer the new coordinates to work array
    2549          83 :             IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
    2550             :                ! everything belongs to the same atom
    2551         316 :                DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
    2552        1066 :                   helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
    2553             :                END DO
    2554             :             ELSE
    2555             :                ! is distributed among two atoms
    2556             :                ! transfer to atom not containing the tail bead
    2557          50 :                DO kbead = helium%worm_bead_idx_work, helium%beads
    2558         149 :                   helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
    2559             :                END DO
    2560             :                ! transfer to atom containing the tail bead
    2561          55 :                DO kbead = 1, helium%worm_bead_idx - 1
    2562         169 :                   helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
    2563             :                END DO
    2564             :             END IF
    2565         332 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2566          83 :             helium%worm_bead_idx_work = helium%worm_bead_idx
    2567          83 :             helium%worm_atom_idx_work = helium%worm_atom_idx
    2568          83 :             ac = 0
    2569         145 :             RETURN
    2570             :          END IF
    2571             :       END IF
    2572             : 
    2573             :       ! for now accepted
    2574             :       ! rejection of the whole move if any centroid is farther away
    2575             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    2576             :       ! AND ist not moving towards the center
    2577        1983 :       IF (.NOT. helium%periodic) THEN
    2578        1969 :          IF (helium%solute_present) THEN
    2579        7876 :             new_com(:) = helium%center(:)
    2580        1969 :             old_com(:) = new_com(:)
    2581             :          ELSE
    2582           0 :             new_com(:) = 0.0_dp
    2583           0 :             DO ia = 1, helium%atoms
    2584           0 :                DO ib = 1, helium%beads
    2585           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    2586             :                END DO
    2587             :             END DO
    2588           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2589             :             ! also compute the old center of mass (optimization potential)
    2590           0 :             old_com(:) = 0.0_dp
    2591           0 :             DO ia = 1, helium%atoms
    2592           0 :                DO ib = 1, helium%beads
    2593           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    2594             :                END DO
    2595             :             END DO
    2596           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2597             :          END IF
    2598        1969 :          should_reject = .FALSE.
    2599       58726 :          atom: DO ia = 1, helium%atoms
    2600       56819 :             dr(:) = 0.0_dp
    2601      965923 :             DO ib = 1, helium%beads
    2602     3693235 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    2603             :             END DO
    2604      227276 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    2605      227276 :             rtmp = DOT_PRODUCT(dr, dr)
    2606       58726 :             IF (rtmp >= helium%droplet_radius**2) THEN
    2607          62 :                dro(:) = 0.0_dp
    2608        1054 :                DO ib = 1, helium%beads
    2609        4030 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    2610             :                END DO
    2611         248 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    2612         248 :                rtmpo = DOT_PRODUCT(dro, dro)
    2613             :                ! only reject if it moves away from COM outside the droplet radius
    2614          62 :                IF (rtmpo < rtmp) THEN
    2615             :                   ! found - this one does not fit within R from the center
    2616             :                   should_reject = .TRUE.
    2617             :                   EXIT atom
    2618             :                END IF
    2619             :             END IF
    2620             :          END DO atom
    2621        1969 :          IF (should_reject) THEN
    2622             :             ! restore original coordinates
    2623             :             ! write back only changed atoms
    2624             :             ! write back only changed atoms
    2625             :             ! transfer the new coordinates to work array
    2626          62 :             IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
    2627             :                ! everything belongs to the same atom
    2628         282 :                DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
    2629         948 :                   helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
    2630             :                END DO
    2631             :             ELSE
    2632             :                ! is distributed among two atoms
    2633             :                ! transfer to atom not containing the tail bead
    2634           7 :                DO kbead = helium%worm_bead_idx_work, helium%beads
    2635          22 :                   helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
    2636             :                END DO
    2637             :                ! transfer to atom containing the tail bead
    2638           2 :                DO kbead = 1, helium%worm_bead_idx - 1
    2639           2 :                   helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
    2640             :                END DO
    2641             :             END IF
    2642         248 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2643          62 :             helium%worm_bead_idx_work = helium%worm_bead_idx
    2644          62 :             helium%worm_atom_idx_work = helium%worm_atom_idx
    2645          62 :             ac = 0
    2646          62 :             RETURN
    2647             :          END IF
    2648             :       END IF
    2649             : 
    2650             :       ! finally accepted
    2651        1921 :       ac = 1
    2652             :       ! write changed coordinates to position array
    2653             :       ! write back only changed atoms
    2654        1921 :       IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
    2655             :          ! everything belongs to the same atom
    2656        6599 :          DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
    2657       21917 :             helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
    2658             :          END DO
    2659             :       ELSE
    2660             :          ! is distributed among two atoms
    2661             :          ! transfer to atom not containing the tail bead
    2662        1465 :          DO kbead = helium%worm_bead_idx_work, helium%beads
    2663        4576 :             helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
    2664             :          END DO
    2665             :          ! transfer to atom containing the tail bead
    2666        1035 :          DO kbead = 1, helium%worm_bead_idx - 1
    2667        2856 :             helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
    2668             :          END DO
    2669             :       END IF
    2670        7684 :       helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
    2671        1921 :       helium%worm_bead_idx = helium%worm_bead_idx_work
    2672        1921 :       helium%worm_atom_idx = helium%worm_atom_idx_work
    2673             : 
    2674             :    END SUBROUTINE worm_crawl_move_backward
    2675             : 
    2676             : ! **************************************************************************************************
    2677             : !> \brief Free particle density matrix
    2678             : !> \param helium ...
    2679             : !> \param startpos ...
    2680             : !> \param endpos ...
    2681             : !> \param l ...
    2682             : !> \return ...
    2683             : !> \author fuhl
    2684             : ! **************************************************************************************************
    2685      596992 :    REAL(KIND=dp) FUNCTION free_density_matrix(helium, startpos, endpos, l) RESULT(rho)
    2686             : 
    2687             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2688             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: startpos, endpos
    2689             :       INTEGER, INTENT(IN)                                :: l
    2690             : 
    2691             :       REAL(KIND=dp)                                      :: distsq, prefac
    2692             :       REAL(KIND=dp), DIMENSION(3)                        :: dvec
    2693             : 
    2694     2387968 :       dvec(:) = startpos(:) - endpos(:)
    2695      596992 :       CALL helium_pbc(helium, dvec)
    2696     2387968 :       distsq = DOT_PRODUCT(dvec, dvec)
    2697             : 
    2698      596992 :       prefac = 0.5_dp/(helium%hb2m*REAL(l, dp)*helium%tau)
    2699      596992 :       rho = -prefac*distsq
    2700      596992 :       rho = EXP(rho)
    2701      596992 :       rho = rho*(SQRT(prefac/pi))**3
    2702             : 
    2703      596992 :    END FUNCTION free_density_matrix
    2704             : 
    2705             : ! **************************************************************************************************
    2706             : !> \brief Swap move in open state to change permutation state
    2707             : !> \param pint_env ...
    2708             : !> \param helium ...
    2709             : !> \param natoms ...
    2710             : !> \param l ...
    2711             : !> \param ac ...
    2712             : !> \author fuhl
    2713             : ! **************************************************************************************************
    2714       10448 :    SUBROUTINE worm_swap_move(pint_env, helium, natoms, l, ac)
    2715             : 
    2716             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2717             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    2718             :       INTEGER, INTENT(IN)                                :: natoms, l
    2719             :       INTEGER, INTENT(OUT)                               :: ac
    2720             : 
    2721             :       INTEGER                                            :: bendatom, bstartatom, endbead, &
    2722             :                                                             excludeatom, fendatom, fstartatom, ia, &
    2723             :                                                             iatom, ib, ibead, jbead, kbead, &
    2724             :                                                             startbead, tmpi
    2725             :       LOGICAL                                            :: should_reject
    2726             :       REAL(KIND=dp)                                      :: backwarddensmatsum, forwarddensmatsum, &
    2727             :                                                             mass, newheadpotential, &
    2728             :                                                             oldheadpotential, rtmp, rtmpo, sdiff, &
    2729             :                                                             snew, sold
    2730             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
    2731       20896 :       REAL(KIND=dp), DIMENSION(3, l)                     :: newsection
    2732       10448 :       REAL(KIND=dp), DIMENSION(natoms)                   :: backwarddensmat, forwarddensmat
    2733             : 
    2734             :       ! first the endbead of the reconstruction is needed
    2735       10448 :       startbead = helium%worm_bead_idx
    2736       10448 :       endbead = helium%worm_bead_idx + l + 1
    2737       10448 :       fstartatom = helium%worm_atom_idx
    2738       10448 :       excludeatom = fstartatom
    2739             :       ! compute the atomwise probabilities to be the worms swap partner
    2740             :       ! Check if the imaginary time section belongs to two atoms
    2741       10448 :       IF (endbead > helium%beads) THEN
    2742        2777 :          endbead = endbead - helium%beads
    2743             :          ! exclude atom is the one not to connect to because it will result in an unnatural state
    2744        2777 :          excludeatom = helium%permutation(excludeatom)
    2745             :       END IF
    2746      155451 :       DO iatom = 1, excludeatom - 1
    2747             :          forwarddensmat(iatom) = free_density_matrix(helium, helium%worm_xtra_bead(:), &
    2748      155451 :                                                      helium%pos(:, iatom, endbead), l + 1)
    2749             :       END DO
    2750       10448 :       forwarddensmat(excludeatom) = 0.0_dp
    2751      163941 :       DO iatom = excludeatom + 1, natoms
    2752             :          forwarddensmat(iatom) = free_density_matrix(helium, helium%worm_xtra_bead(:), &
    2753      163941 :                                                      helium%pos(:, iatom, endbead), l + 1)
    2754             :       END DO
    2755      319392 :       forwarddensmatsum = SUM(forwarddensmat)
    2756       10448 :       IF (forwarddensmatsum == 0.0_dp) THEN
    2757           0 :          ac = 0
    2758           0 :          RETURN
    2759             :       END IF
    2760             : 
    2761             :       ! Select an atom with its corresponding probability
    2762       10448 :       rtmp = helium%rng_stream_uniform%next()*forwarddensmatsum
    2763       10448 :       fendatom = 1
    2764      167857 :       DO WHILE (rtmp >= forwarddensmat(fendatom))
    2765      157409 :          rtmp = rtmp - forwarddensmat(fendatom)
    2766      157409 :          fendatom = fendatom + 1
    2767             :       END DO
    2768             :       ! just for numerical safety
    2769       10448 :       fendatom = MIN(fendatom, natoms)
    2770       10448 :       IF (fendatom == excludeatom) THEN
    2771           0 :          ac = 0
    2772           0 :          RETURN
    2773             :       END IF
    2774             : 
    2775             :       ! ensure detailed balance
    2776       10448 :       IF (endbead > startbead) THEN
    2777        7671 :          bstartatom = fendatom
    2778             :       ELSE
    2779        2777 :          bstartatom = helium%iperm(fendatom)
    2780             :       END IF
    2781             :       bendatom = fendatom
    2782             : 
    2783      155451 :       DO iatom = 1, excludeatom - 1
    2784             :          backwarddensmat(iatom) = free_density_matrix(helium, &
    2785             :                                                       helium%pos(:, bstartatom, startbead), &
    2786      155451 :                                                       helium%pos(:, iatom, endbead), l + 1)
    2787             :       END DO
    2788       10448 :       backwarddensmat(excludeatom) = 0.0_dp
    2789      163941 :       DO iatom = excludeatom + 1, natoms
    2790             :          backwarddensmat(iatom) = free_density_matrix(helium, &
    2791             :                                                       helium%pos(:, bstartatom, startbead), &
    2792      163941 :                                                       helium%pos(:, iatom, endbead), l + 1)
    2793             :       END DO
    2794             : 
    2795      319392 :       backwarddensmatsum = SUM(backwarddensmat)
    2796             : 
    2797       10448 :       mass = helium%he_mass_au
    2798             : 
    2799             :       !compute action before the move
    2800             :       sold = worm_path_action(helium, helium%pos, &
    2801       10448 :                               bstartatom, startbead, fendatom, endbead)
    2802             : 
    2803       10448 :       IF (helium%solute_present) THEN
    2804             :          ! no special head treatment needed, as it will change due to swapping
    2805             :          ! the worm gap and due to primitive coupling no cross bead terms are considered
    2806             :          sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
    2807       10268 :                                               bstartatom, startbead, fendatom, endbead)
    2808             :          ! compute potential of old and new head here (only once, as later only a rescaling is necessary)
    2809             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2810             :                                      fstartatom, startbead, helium%worm_xtra_bead, &
    2811       10268 :                                      energy=oldheadpotential)
    2812       10268 :          oldheadpotential = oldheadpotential*helium%tau
    2813             : 
    2814             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    2815             :                                      bstartatom, startbead, helium%pos(:, bstartatom, startbead), &
    2816       10268 :                                      energy=newheadpotential)
    2817       10268 :          newheadpotential = newheadpotential*helium%tau
    2818             : 
    2819       10268 :          sold = sold + 0.5_dp*oldheadpotential + newheadpotential
    2820             :       END IF
    2821             : 
    2822             :       ! construct a new path connecting the start and endbead
    2823             :       ! need to be the old coordinates due to reordering of the work coordinates
    2824             :       CALL path_construct(helium, &
    2825             :                           helium%worm_xtra_bead(:), &
    2826             :                           helium%pos(:, fendatom, endbead), l, &
    2827       10448 :                           newsection)
    2828             : 
    2829             :       !write new path segment to work array
    2830             :       !first the part that is guaranteed to fit on the coorinates of startatom of the
    2831       10448 :       jbead = 1
    2832       10448 :       IF (startbead < endbead) THEN
    2833             :          ! everything belongs to the same atom
    2834       33777 :          DO kbead = startbead + 1, endbead - 1
    2835      104424 :             helium%work(:, fstartatom, kbead) = newsection(:, jbead)
    2836       33777 :             jbead = jbead + 1
    2837             :          END DO
    2838             :       ELSE
    2839             :          ! is distributed among two atoms
    2840             :          ! transfer to atom not containing the head bead
    2841        7977 :          DO kbead = startbead + 1, helium%beads
    2842       20800 :             helium%work(:, fstartatom, kbead) = newsection(:, jbead)
    2843        7977 :             jbead = jbead + 1
    2844             :          END DO
    2845             :          ! rest to the second atom
    2846        8017 :          DO ibead = 1, endbead - 1
    2847       20960 :             helium%work(:, fendatom, ibead) = newsection(:, jbead)
    2848        8017 :             jbead = jbead + 1
    2849             :          END DO
    2850             :       END IF
    2851             : 
    2852             :       !exchange coordinates head coordinate first
    2853       41792 :       helium%work(:, helium%worm_atom_idx, helium%worm_bead_idx) = helium%worm_xtra_bead(:)
    2854       41792 :       helium%worm_xtra_bead_work(:) = helium%pos(:, bstartatom, startbead)
    2855             : 
    2856             :       ! move coordinates from former worm atom tail to new worm atom tail
    2857      102032 :       DO ib = startbead, helium%beads
    2858      376784 :          helium%work(:, bstartatom, ib) = helium%pos(:, fstartatom, ib)
    2859             :       END DO
    2860             :       ! need to copy the rest of pselatom to former worm atom
    2861       10448 :       IF (endbead > startbead) THEN
    2862       57501 :          DO ib = endbead, helium%beads
    2863      206991 :             helium%work(:, fstartatom, ib) = helium%pos(:, bstartatom, ib)
    2864             :          END DO
    2865             :       END IF
    2866             : 
    2867             :       !update permtable
    2868       10448 :       tmpi = helium%permutation(bstartatom)
    2869       10448 :       helium%permutation(bstartatom) = helium%permutation(fstartatom)
    2870       10448 :       helium%permutation(fstartatom) = tmpi
    2871             :       !update inverse permtable
    2872      319392 :       DO ib = 1, SIZE(helium%permutation)
    2873      319392 :          helium%iperm(helium%permutation(ib)) = ib
    2874             :       END DO
    2875       10448 :       helium%worm_atom_idx_work = bstartatom
    2876             :       ! action after the move
    2877             : 
    2878       10448 :       IF (endbead > startbead) THEN
    2879             :          snew = worm_path_action(helium, helium%work, &
    2880        7671 :                                  fstartatom, startbead, fstartatom, endbead)
    2881        7671 :          IF (helium%solute_present) THEN
    2882             :             ! no special head treatment needed, because a swap can't go over
    2883             :             ! the worm gap and due to primitive coupling no cross bead terms are considered
    2884             :             snew = snew + worm_path_inter_action(pint_env, helium, helium%work, &
    2885        7577 :                                                  fstartatom, startbead, fstartatom, endbead)
    2886             : 
    2887             :             ! add the previously computed old and new head actions
    2888        7577 :             snew = snew + oldheadpotential + 0.5_dp*newheadpotential
    2889             :          END IF
    2890             :       ELSE
    2891             :          snew = worm_path_action(helium, helium%work, &
    2892        2777 :                                  fstartatom, startbead, helium%permutation(fstartatom), endbead)
    2893        2777 :          IF (helium%solute_present) THEN
    2894             :             ! no special head treatment needed, because a swap can't go over
    2895             :             ! the worm gap and due to primitive coupling no cross bead terms are considered
    2896             :             snew = snew + worm_path_inter_action(pint_env, helium, helium%work, &
    2897        2691 :                                                  fstartatom, startbead, helium%permutation(fstartatom), endbead)
    2898             : 
    2899             :             ! add the previously computed old and new head actions
    2900        2691 :             snew = snew + oldheadpotential + 0.5_dp*newheadpotential
    2901             :          END IF
    2902             :       END IF
    2903             : 
    2904             :       ! Metropolis:
    2905       10448 :       sdiff = sold - snew
    2906       10448 :       sdiff = sdiff + LOG(forwarddensmatsum/backwarddensmatsum)
    2907       10448 :       IF (sdiff < 0) THEN
    2908       10424 :          should_reject = .FALSE.
    2909       10424 :          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
    2910             :             should_reject = .TRUE.
    2911             :          ELSE
    2912        4792 :             rtmp = helium%rng_stream_uniform%next()
    2913        4792 :             IF (EXP(sdiff) < rtmp) THEN
    2914             :                should_reject = .TRUE.
    2915             :             END IF
    2916             :          END IF
    2917             :          IF (should_reject) THEN
    2918             :             ! rejected !
    2919             :             ! write back only changed atoms
    2920      101550 :             DO kbead = startbead, helium%beads
    2921      374970 :                helium%work(:, bstartatom, kbead) = helium%pos(:, bstartatom, kbead)
    2922             :             END DO
    2923      101550 :             DO kbead = startbead, helium%beads
    2924      374970 :                helium%work(:, fstartatom, kbead) = helium%pos(:, fstartatom, kbead)
    2925             :             END DO
    2926       41640 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    2927       10410 :             helium%worm_atom_idx_work = helium%worm_atom_idx
    2928       10410 :             IF (startbead > endbead) THEN
    2929       10794 :                DO kbead = 1, endbead
    2930       34845 :                   helium%work(:, fendatom, kbead) = helium%pos(:, fendatom, kbead)
    2931             :                END DO
    2932             :             END IF
    2933             :             ! reset permtable
    2934       10410 :             tmpi = helium%permutation(bstartatom)
    2935       10410 :             helium%permutation(bstartatom) = helium%permutation(fstartatom)
    2936       10410 :             helium%permutation(fstartatom) = tmpi
    2937             :             !update inverse permtable
    2938      318138 :             DO ib = 1, SIZE(helium%permutation)
    2939      318138 :                helium%iperm(helium%permutation(ib)) = ib
    2940             :             END DO
    2941       10410 :             ac = 0
    2942       10410 :             RETURN
    2943             :          END IF
    2944             :       END IF
    2945             : 
    2946             :       ! for now accepted
    2947             :       ! rejection of the whole move if any centroid is farther away
    2948             :       ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
    2949             :       ! AND ist not moving towards the center
    2950          38 :       IF (.NOT. helium%periodic) THEN
    2951          38 :          IF (helium%solute_present) THEN
    2952         152 :             new_com(:) = helium%center(:)
    2953          38 :             old_com(:) = new_com(:)
    2954             :          ELSE
    2955           0 :             new_com(:) = 0.0_dp
    2956           0 :             DO ia = 1, helium%atoms
    2957           0 :                DO ib = 1, helium%beads
    2958           0 :                   new_com(:) = new_com(:) + helium%work(:, ia, ib)
    2959             :                END DO
    2960             :             END DO
    2961           0 :             new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2962             :             ! also compute the old center of mass (optimization potential)
    2963           0 :             old_com(:) = 0.0_dp
    2964           0 :             DO ia = 1, helium%atoms
    2965           0 :                DO ib = 1, helium%beads
    2966           0 :                   old_com(:) = old_com(:) + helium%pos(:, ia, ib)
    2967             :                END DO
    2968             :             END DO
    2969           0 :             old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
    2970             :          END IF
    2971          38 :          should_reject = .FALSE.
    2972        1074 :          atom: DO ia = 1, helium%atoms
    2973        1066 :             dr(:) = 0.0_dp
    2974       18122 :             DO ib = 1, helium%beads
    2975       69290 :                dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
    2976             :             END DO
    2977        4264 :             dr(:) = dr(:)/REAL(helium%beads, dp)
    2978        4264 :             rtmp = DOT_PRODUCT(dr, dr)
    2979        1074 :             IF (rtmp >= helium%droplet_radius**2) THEN
    2980          30 :                dro(:) = 0.0_dp
    2981         510 :                DO ib = 1, helium%beads
    2982        1950 :                   dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
    2983             :                END DO
    2984         120 :                dro(:) = dro(:)/REAL(helium%beads, dp)
    2985         120 :                rtmpo = DOT_PRODUCT(dro, dro)
    2986             :                ! only reject if it moves away from COM outside the droplet radius
    2987          30 :                IF (rtmpo < rtmp) THEN
    2988             :                   ! found - this one does not fit within R from the center
    2989             :                   should_reject = .TRUE.
    2990             :                   EXIT atom
    2991             :                END IF
    2992             :             END IF
    2993             :          END DO atom
    2994          38 :          IF (should_reject) THEN
    2995             :             ! write back only changed atoms
    2996         380 :             DO kbead = startbead, helium%beads
    2997        1430 :                helium%work(:, bstartatom, kbead) = helium%pos(:, bstartatom, kbead)
    2998             :             END DO
    2999         380 :             DO kbead = startbead, helium%beads
    3000        1430 :                helium%work(:, fstartatom, kbead) = helium%pos(:, fstartatom, kbead)
    3001             :             END DO
    3002         120 :             helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
    3003          30 :             helium%worm_atom_idx_work = helium%worm_atom_idx
    3004          30 :             IF (startbead > endbead) THEN
    3005           0 :                DO kbead = 1, endbead
    3006           0 :                   helium%work(:, fendatom, kbead) = helium%pos(:, fendatom, kbead)
    3007             :                END DO
    3008             :             END IF
    3009             : 
    3010             :             ! reset permtable
    3011          30 :             tmpi = helium%permutation(bstartatom)
    3012          30 :             helium%permutation(bstartatom) = helium%permutation(fstartatom)
    3013          30 :             helium%permutation(fstartatom) = tmpi
    3014             :             !update inverse permtable
    3015         990 :             DO ib = 1, SIZE(helium%permutation)
    3016         990 :                helium%iperm(helium%permutation(ib)) = ib
    3017             :             END DO
    3018          30 :             ac = 0
    3019          30 :             RETURN
    3020             :          END IF
    3021             :       END IF
    3022             : 
    3023             :       ! finally accepted
    3024           8 :       ac = 1
    3025         102 :       DO kbead = startbead, helium%beads
    3026         384 :          helium%pos(:, bstartatom, kbead) = helium%work(:, bstartatom, kbead)
    3027             :       END DO
    3028         102 :       DO kbead = startbead, helium%beads
    3029         384 :          helium%pos(:, fstartatom, kbead) = helium%work(:, fstartatom, kbead)
    3030             :       END DO
    3031          32 :       helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
    3032           8 :       helium%worm_atom_idx = helium%worm_atom_idx_work
    3033           8 :       IF (startbead > endbead) THEN
    3034           0 :          DO kbead = 1, endbead
    3035           0 :             helium%pos(:, fendatom, kbead) = helium%work(:, fendatom, kbead)
    3036             :          END DO
    3037             :       END IF
    3038             : 
    3039             :    END SUBROUTINE worm_swap_move
    3040             : 
    3041             : ! **************************************************************************************************
    3042             : !> \brief Action along path
    3043             : !> \param helium ...
    3044             : !> \param pos ...
    3045             : !> \param startatom ...
    3046             : !> \param startbead ...
    3047             : !> \param endatom ...
    3048             : !> \param endbead ...
    3049             : !> \return ...
    3050             : !> \author fuhl
    3051             : ! **************************************************************************************************
    3052       42167 :    REAL(KIND=dp) FUNCTION worm_path_action(helium, pos, &
    3053             :                                            startatom, startbead, endatom, endbead) RESULT(partaction)
    3054             : 
    3055             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    3056             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
    3057             :          POINTER                                         :: pos
    3058             :       INTEGER, INTENT(IN)                                :: startatom, startbead, endatom, endbead
    3059             : 
    3060             :       INTEGER                                            :: iatom, ibead
    3061       42167 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work2, work3
    3062       42167 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
    3063             : 
    3064      295169 :       ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
    3065       42167 :       partaction = 0.0_dp
    3066             :       ! do action in two ways, depending
    3067             :       ! if coordinates are not wrapping
    3068       42167 :       IF (startbead < endbead) THEN
    3069             :          ! helium pair action
    3070             :          ! every atom, with the one (or two) who got a resampling
    3071      945045 :          DO iatom = 1, helium%atoms
    3072             :             ! avoid self interaction
    3073      913844 :             IF (iatom == startatom) CYCLE
    3074             :             ! first the section up to the worm gap
    3075             :             ! two less, because we need to work on the worm intersection separately
    3076     5452076 :             DO ibead = startbead, endbead
    3077    19160375 :                work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3078             :             END DO
    3079      945045 :             partaction = partaction + helium_eval_chain(helium, work, endbead - startbead + 1, work2, work3)
    3080             :          END DO
    3081             :       ELSE !(startbead > endbead)
    3082             :          ! helium pair action
    3083             :          ! every atom, with the one (or two) who got a resampling
    3084      329106 :          DO iatom = 1, helium%atoms
    3085             :             ! avoid self interaction
    3086      318140 :             IF (iatom == startatom) CYCLE
    3087             :             ! first the section up to the worm gap
    3088             :             ! two less, because we need to work on the worm intersection separately
    3089     1171913 :             DO ibead = startbead, helium%beads
    3090     3766130 :                work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3091             :             END DO
    3092             :             ! wrapping bond
    3093     1228696 :             work(:, helium%beads + 2 - startbead) = pos(:, helium%permutation(iatom), 1) - pos(:, endatom, 1)
    3094      329106 :             partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 2, work2, work3)
    3095             :          END DO
    3096             : 
    3097             :          ! ending atom
    3098      329106 :          DO iatom = 1, helium%atoms
    3099             :             ! avoid self interaction
    3100      318140 :             IF (iatom == endatom) CYCLE
    3101             :             !from first to endbead
    3102      318140 :             IF (endbead > 1) THEN
    3103     1021303 :                DO ibead = 1, endbead
    3104     3377137 :                   work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
    3105             :                END DO
    3106      236025 :                partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
    3107             :             END IF
    3108             :          END DO
    3109             : 
    3110             :       END IF
    3111       42167 :       DEALLOCATE (work, work2, work3)
    3112             : 
    3113       42167 :    END FUNCTION worm_path_action
    3114             : 
    3115             : ! **************************************************************************************************
    3116             : !> \brief Corrected path action for worm
    3117             : !> \param helium ...
    3118             : !> \param pos ...
    3119             : !> \param startatom ...
    3120             : !> \param startbead ...
    3121             : !> \param endatom ...
    3122             : !> \param endbead ...
    3123             : !> \param xtrapos ...
    3124             : !> \param worm_atom_idx ...
    3125             : !> \param worm_bead_idx ...
    3126             : !> \return ...
    3127             : !> \author fuhl
    3128             : ! **************************************************************************************************
    3129        9479 :    REAL(KIND=dp) FUNCTION worm_path_action_worm_corrected(helium, pos, &
    3130             :                                    startatom, startbead, endatom, endbead, xtrapos, worm_atom_idx, worm_bead_idx) RESULT(partaction)
    3131             : 
    3132             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    3133             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
    3134             :          POINTER                                         :: pos
    3135             :       INTEGER, INTENT(IN)                                :: startatom, startbead, endatom, endbead
    3136             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: xtrapos
    3137             :       INTEGER, INTENT(IN)                                :: worm_atom_idx, worm_bead_idx
    3138             : 
    3139             :       INTEGER                                            :: iatom, ibead
    3140        9479 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work2, work3
    3141        9479 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
    3142             :       REAL(KIND=dp), DIMENSION(3)                        :: r, rp
    3143             : 
    3144       66353 :       ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
    3145        9479 :       partaction = 0.0_dp
    3146             :       ! do action in two ways, depending
    3147             :       ! if coordinates are not wrapping
    3148        9479 :       IF (startbead < endbead) THEN
    3149             :          ! helium pair action
    3150             :          ! every atom, with the one (or two) who got a resampling
    3151      222417 :          DO iatom = 1, helium%atoms
    3152             :             ! avoid self interaction
    3153      215124 :             IF (iatom == startatom) CYCLE
    3154             :             ! first the section up to the worm gap
    3155             :             ! two less, because we need to work on the worm intersection separately
    3156      207831 :             IF (worm_bead_idx - startbead > 1) THEN
    3157      906831 :                DO ibead = startbead, worm_bead_idx - 1
    3158     3013527 :                   work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3159             :                END DO
    3160      204599 :                partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - startbead, work2, work3)
    3161             :             END IF
    3162             : 
    3163      215124 :             IF (endbead > worm_bead_idx) THEN
    3164       44200 :                DO ibead = worm_bead_idx, endbead
    3165      146452 :                   work(:, ibead + 1 - worm_bead_idx) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3166             :                END DO
    3167       10116 :                partaction = partaction + helium_eval_chain(helium, work, endbead - worm_bead_idx + 1, work2, work3)
    3168             :             END IF
    3169             :          END DO
    3170             : 
    3171        7293 :          IF (worm_atom_idx /= startatom) THEN
    3172       13572 :             DO iatom = 1, helium%atoms
    3173       13136 :                IF (iatom == startatom) CYCLE
    3174       12700 :                IF (iatom == worm_atom_idx) CYCLE
    3175       49056 :                r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
    3176       49056 :                rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, startatom, worm_bead_idx)
    3177       13572 :                partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3178             :             END DO
    3179             :             ! add pair action with worm
    3180        1744 :             r(:) = pos(:, startatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
    3181        1744 :             rp(:) = pos(:, startatom, worm_bead_idx) - xtrapos(:)
    3182         436 :             partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3183             :          ELSE
    3184             :             ! worm intersection
    3185      208845 :             DO iatom = 1, helium%atoms
    3186      201988 :                IF (iatom == startatom) CYCLE
    3187      780524 :                r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
    3188      780524 :                rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
    3189      208845 :                partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3190             :             END DO
    3191             :          END IF
    3192             :       ELSE !(startbead > endbead)
    3193             :          ! section wraps around the atom
    3194             :          ! two cases: worm gap is before or after wrap
    3195        2186 :          IF (worm_bead_idx > startbead) THEN
    3196             :             ! action terms up to worm gap on starting atom
    3197        1656 :             DO iatom = 1, helium%atoms
    3198             :                ! avoid self interaction
    3199        1600 :                IF (iatom == startatom) CYCLE
    3200             :                ! first the section up to the worm gap
    3201             :                ! two less, because we need to work on the worm intersection separately
    3202        1544 :                IF (worm_bead_idx - startbead > 1) THEN
    3203        2450 :                   DO ibead = startbead, worm_bead_idx - 1
    3204        7598 :                      work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3205             :                   END DO
    3206         734 :                   partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - startbead, work2, work3)
    3207             :                END IF
    3208             : 
    3209             :                ! up to the wrapping border
    3210        4374 :                DO ibead = worm_bead_idx, helium%beads
    3211       12864 :                   work(:, ibead + 1 - worm_bead_idx) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3212             :                END DO
    3213             :                ! wrapping bond
    3214             :                work(:, helium%beads - worm_bead_idx + 2) = pos(:, helium%permutation(iatom), 1) - &
    3215        6176 :                                                            pos(:, helium%permutation(startatom), 1)
    3216        1656 :                partaction = partaction + helium_eval_chain(helium, work, helium%beads - worm_bead_idx + 2, work2, work3)
    3217             : 
    3218             :             END DO
    3219             : 
    3220             :             ! ending atom
    3221        1656 :             DO iatom = 1, helium%atoms
    3222             :                ! avoid self interaction
    3223        1600 :                IF (iatom == endatom) CYCLE
    3224             :                !from first to endbead
    3225        1600 :                IF (endbead > 1) THEN
    3226        3202 :                   DO ibead = 1, endbead
    3227       10006 :                      work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
    3228             :                   END DO
    3229         934 :                   partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
    3230             :                END IF
    3231             :             END DO
    3232             : 
    3233          56 :             IF (worm_atom_idx /= startatom) THEN
    3234        1656 :                DO iatom = 1, helium%atoms
    3235        1600 :                   IF (iatom == startatom) CYCLE
    3236        1544 :                   IF (iatom == worm_atom_idx) CYCLE
    3237        5952 :                   r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
    3238        5952 :                   rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, startatom, worm_bead_idx)
    3239        1656 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3240             :                END DO
    3241             :                ! add pair action with worm
    3242         224 :                r(:) = pos(:, startatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
    3243         224 :                rp(:) = pos(:, startatom, worm_bead_idx) - xtrapos(:)
    3244          56 :                partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3245             :             ELSE
    3246             :                ! worm intersection
    3247           0 :                DO iatom = 1, helium%atoms
    3248           0 :                   IF (iatom == startatom) CYCLE
    3249           0 :                   r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
    3250           0 :                   rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
    3251           0 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3252             :                END DO
    3253             :             END IF
    3254             :          ELSE !(worm_bead_idx < endbead)
    3255        2130 :             IF (worm_bead_idx /= 1) THEN
    3256       44553 :                DO iatom = 1, helium%atoms
    3257             :                   ! avoid self interaction
    3258       43072 :                   IF (iatom == startatom) CYCLE
    3259             :                   ! first the section up to the end of the atom
    3260             :                   ! one less, because we need to work on the wrapping separately
    3261      124787 :                   DO ibead = startbead, helium%beads
    3262      374375 :                      work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3263             :                   END DO
    3264             :                   ! wrapping bond
    3265      166364 :              work(:, helium%beads - startbead + 2) = pos(:, helium%permutation(iatom), 1) - pos(:, helium%permutation(startatom), 1)
    3266       44553 :                   partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 2, work2, work3)
    3267             :                END DO
    3268             : 
    3269             :                ! ending atom
    3270       44553 :                DO iatom = 1, helium%atoms
    3271             :                   ! avoid self interaction
    3272       43072 :                   IF (iatom == endatom) CYCLE
    3273             :                   !from first to two before the worm gap
    3274       41591 :                   IF (worm_bead_idx > 2) THEN
    3275       90118 :                      DO ibead = 1, worm_bead_idx - 1
    3276      285934 :                         work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
    3277             :                      END DO
    3278       24846 :                      partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - 1, work2, work3)
    3279             :                   END IF
    3280             : 
    3281       43072 :                   IF (endbead > worm_bead_idx) THEN
    3282        5484 :                      DO ibead = worm_bead_idx, endbead
    3283       17358 :                         work(:, ibead - worm_bead_idx + 1) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
    3284             :                      END DO
    3285        1526 :                      partaction = partaction + helium_eval_chain(helium, work, endbead - worm_bead_idx + 1, work2, work3)
    3286             :                   END IF
    3287             :                END DO
    3288             : 
    3289        1481 :                IF (worm_atom_idx /= endatom) THEN
    3290        2016 :                   DO iatom = 1, helium%atoms
    3291        1952 :                      IF (iatom == endatom) CYCLE
    3292        1888 :                      IF (iatom == worm_atom_idx) CYCLE
    3293        7296 :                      r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, endatom, worm_bead_idx - 1)
    3294        7296 :                      rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, endatom, worm_bead_idx)
    3295        2016 :                      partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3296             :                   END DO
    3297             :                   ! add pair action with worm
    3298         256 :                   r(:) = pos(:, endatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
    3299         256 :                   rp(:) = pos(:, endatom, worm_bead_idx) - xtrapos(:)
    3300          64 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3301             :                ELSE
    3302             :                   ! worm intersection
    3303       42537 :                   DO iatom = 1, helium%atoms
    3304       41120 :                      IF (iatom == endatom) CYCLE
    3305      158812 :                      r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, endatom, worm_bead_idx - 1)
    3306      158812 :                      rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
    3307       42537 :                      partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3308             :                   END DO
    3309             :                END IF
    3310             :             ELSE !(worm_bead_idx == 1)
    3311             :                ! special treatment if first bead is worm bead
    3312       19917 :                DO iatom = 1, helium%atoms
    3313             :                   ! avoid self interaction
    3314       19268 :                   IF (iatom == startatom) CYCLE
    3315             :                   ! first the section up to the end of the atom
    3316             :                   ! one less, because we need to work on the wrapping separately
    3317       19268 :                   IF (helium%beads > startbead) THEN
    3318       79426 :                      DO ibead = startbead, helium%beads
    3319      263335 :                         work(:, ibead - startbead + 1) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
    3320             :                      END DO
    3321       18123 :                      partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 1, work2, work3)
    3322             :                   END IF
    3323             :                END DO
    3324             : 
    3325             :                ! ending atom
    3326       19917 :                DO iatom = 1, helium%atoms
    3327       19917 :                   IF (endbead > 1) THEN
    3328        6144 :                      DO ibead = 1, endbead
    3329       20736 :                         work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
    3330             :                      END DO
    3331        1280 :                      partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
    3332             :                   END IF
    3333             :                END DO
    3334             : 
    3335         649 :                IF (worm_atom_idx /= endatom) THEN
    3336        1626 :                   DO iatom = 1, helium%atoms
    3337        1576 :                      IF (iatom == endatom) CYCLE
    3338        1526 :                      IF (iatom == worm_atom_idx) CYCLE
    3339        5904 :                      r(:) = pos(:, helium%iperm(iatom), helium%beads) - pos(:, startatom, helium%beads)
    3340        5904 :                      rp(:) = pos(:, iatom, 1) - pos(:, endatom, 1)
    3341        1626 :                      partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3342             :                   END DO
    3343             :                   ! add pair action with worm
    3344         200 :                   r(:) = pos(:, startatom, helium%beads) - pos(:, helium%iperm(worm_atom_idx), helium%beads)
    3345         200 :                   rp(:) = pos(:, endatom, 1) - xtrapos(:)
    3346          50 :                   partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3347             :                ELSE
    3348             :                   ! worm intersection
    3349       18291 :                   DO iatom = 1, helium%atoms
    3350       17692 :                      IF (iatom == endatom) CYCLE
    3351       68372 :                      r(:) = pos(:, helium%iperm(iatom), helium%beads) - pos(:, startatom, helium%beads)
    3352       68372 :                      rp(:) = pos(:, iatom, 1) - xtrapos(:)
    3353       18291 :                      partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
    3354             :                   END DO
    3355             :                END IF
    3356             :             END IF
    3357             :          END IF
    3358             :       END IF
    3359        9479 :       DEALLOCATE (work, work2, work3)
    3360             : 
    3361        9479 :    END FUNCTION worm_path_action_worm_corrected
    3362             : 
    3363             : ! **************************************************************************************************
    3364             : !> \brief Path interaction
    3365             : !> \param pint_env ...
    3366             : !> \param helium ...
    3367             : !> \param pos ...
    3368             : !> \param startatom ...
    3369             : !> \param startbead ...
    3370             : !> \param endatom ...
    3371             : !> \param endbead ...
    3372             : !> \return ...
    3373             : !> \author fuhl
    3374             : ! **************************************************************************************************
    3375       36998 :    REAL(KIND=dp) FUNCTION worm_path_inter_action(pint_env, helium, pos, &
    3376             :                                                  startatom, startbead, endatom, endbead) RESULT(partaction)
    3377             : 
    3378             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    3379             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    3380             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
    3381             :          POINTER                                         :: pos
    3382             :       INTEGER, INTENT(IN)                                :: startatom, startbead, endatom, endbead
    3383             : 
    3384             :       INTEGER                                            :: ibead
    3385             :       REAL(KIND=dp)                                      :: energy
    3386             : 
    3387       36998 :       partaction = 0.0_dp
    3388             : 
    3389             :       ! helium interaction with the solute
    3390             :       ! do action in two ways, depending
    3391             :       ! if coordinates are not wrapping
    3392       36998 :       IF (startbead < endbead) THEN
    3393             : 
    3394             :          ! interaction is only beadwise due to primitive coupling
    3395             :          ! startatom and endatom are the same
    3396      116706 :          DO ibead = startbead + 1, endbead - 1
    3397             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3398       89515 :                                         startatom, ibead, pos(:, startatom, ibead), energy=energy)
    3399      116706 :             partaction = partaction + energy
    3400             :          END DO
    3401             : 
    3402             :       ELSE !(startbead > endbead)
    3403             : 
    3404             :          ! interaction is only beadwise due to primitive coupling
    3405             :          ! startatom and endatom are different
    3406       27909 :          DO ibead = startbead + 1, helium%beads
    3407             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3408       18102 :                                         startatom, ibead, pos(:, startatom, ibead), energy=energy)
    3409       27909 :             partaction = partaction + energy
    3410             :          END DO
    3411             :          ! second atom after imaginary time wrap
    3412       27906 :          DO ibead = 1, endbead - 1
    3413             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3414       18099 :                                         endatom, ibead, pos(:, endatom, ibead), energy=energy)
    3415       27906 :             partaction = partaction + energy
    3416             :          END DO
    3417             :       END IF
    3418             : 
    3419       36998 :       partaction = partaction*helium%tau
    3420             : 
    3421       36998 :    END FUNCTION worm_path_inter_action
    3422             : 
    3423             : ! **************************************************************************************************
    3424             : !> \brief Path interaction for head
    3425             : !> \param pint_env ...
    3426             : !> \param helium ...
    3427             : !> \param pos ...
    3428             : !> \param startatom ...
    3429             : !> \param startbead ...
    3430             : !> \param xtrapos ...
    3431             : !> \param worm_atom_idx ...
    3432             : !> \param worm_bead_idx ...
    3433             : !> \return ...
    3434             : !> \author fuhl
    3435             : ! **************************************************************************************************
    3436        9502 :    REAL(KIND=dp) FUNCTION worm_path_inter_action_head(pint_env, helium, pos, &
    3437             :                                                      startatom, startbead, xtrapos, worm_atom_idx, worm_bead_idx) RESULT(partaction)
    3438             : 
    3439             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    3440             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    3441             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
    3442             :          POINTER                                         :: pos
    3443             :       INTEGER, INTENT(IN)                                :: startatom, startbead
    3444             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: xtrapos
    3445             :       INTEGER, INTENT(IN)                                :: worm_atom_idx, worm_bead_idx
    3446             : 
    3447             :       INTEGER                                            :: ibead
    3448             :       REAL(KIND=dp)                                      :: energy
    3449             : 
    3450        9502 :       partaction = 0.0_dp
    3451             : 
    3452             :       ! helium interaction with the solute
    3453             :       ! if coordinates are not wrapping
    3454        9502 :       IF (startbead < worm_bead_idx) THEN
    3455       24866 :          DO ibead = startbead + 1, worm_bead_idx - 1
    3456             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3457       17602 :                                         startatom, ibead, pos(:, startatom, ibead), energy=energy)
    3458       24866 :             partaction = partaction + energy
    3459             :          END DO
    3460             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    3461        7264 :                                      startatom, ibead, xtrapos, energy=energy)
    3462        7264 :          partaction = partaction + 0.5_dp*energy
    3463             :       ELSE !(startbead > worm_bead_idx)
    3464        5378 :          DO ibead = startbead + 1, helium%beads
    3465             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3466        3140 :                                         startatom, ibead, pos(:, startatom, ibead), energy=energy)
    3467        5378 :             partaction = partaction + energy
    3468             :          END DO
    3469        5360 :          DO ibead = 1, worm_bead_idx - 1
    3470             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3471        3122 :                                         worm_atom_idx, ibead, pos(:, worm_atom_idx, ibead), energy=energy)
    3472        5360 :             partaction = partaction + energy
    3473             :          END DO
    3474             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    3475        2238 :                                      worm_atom_idx, ibead, xtrapos, energy=energy)
    3476        2238 :          partaction = partaction + 0.5_dp*energy
    3477             :       END IF
    3478             : 
    3479        9502 :       partaction = partaction*helium%tau
    3480             : 
    3481        9502 :    END FUNCTION worm_path_inter_action_head
    3482             : 
    3483             : ! **************************************************************************************************
    3484             : !> \brief Path interaction for tail
    3485             : !> \param pint_env ...
    3486             : !> \param helium ...
    3487             : !> \param pos ...
    3488             : !> \param endatom ...
    3489             : !> \param endbead ...
    3490             : !> \param worm_atom_idx ...
    3491             : !> \param worm_bead_idx ...
    3492             : !> \return ...
    3493             : !> \author fuhl
    3494             : ! **************************************************************************************************
    3495        4090 :    REAL(KIND=dp) FUNCTION worm_path_inter_action_tail(pint_env, helium, pos, &
    3496             :                                                       endatom, endbead, worm_atom_idx, worm_bead_idx) RESULT(partaction)
    3497             : 
    3498             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    3499             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    3500             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
    3501             :          POINTER                                         :: pos
    3502             :       INTEGER, INTENT(IN)                                :: endatom, endbead, worm_atom_idx, &
    3503             :                                                             worm_bead_idx
    3504             : 
    3505             :       INTEGER                                            :: ibead
    3506             :       REAL(KIND=dp)                                      :: energy
    3507             : 
    3508        4090 :       partaction = 0.0_dp
    3509             : 
    3510             :       ! helium interaction with the solute
    3511             :       ! if coordinates are not wrapping
    3512        4090 :       IF (worm_bead_idx < endbead) THEN
    3513             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    3514        3299 :                                      worm_atom_idx, worm_bead_idx, pos(:, worm_atom_idx, worm_bead_idx), energy=energy)
    3515        3299 :          partaction = partaction + 0.5_dp*energy
    3516       11260 :          DO ibead = worm_bead_idx + 1, endbead - 1
    3517             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3518        7961 :                                         endatom, ibead, pos(:, endatom, ibead), energy=energy)
    3519       11260 :             partaction = partaction + energy
    3520             :          END DO
    3521             :       ELSE !(worm_bead_idx > endbead)
    3522             :          CALL helium_bead_solute_e_f(pint_env, helium, &
    3523         791 :                                      worm_atom_idx, worm_bead_idx, pos(:, worm_atom_idx, worm_bead_idx), energy=energy)
    3524         791 :          partaction = partaction + 0.5_dp*energy
    3525        2001 :          DO ibead = worm_bead_idx + 1, helium%beads
    3526             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3527        1210 :                                         worm_atom_idx, ibead, pos(:, worm_atom_idx, ibead), energy=energy)
    3528        2001 :             partaction = partaction + energy
    3529             :          END DO
    3530        1818 :          DO ibead = 1, endbead - 1
    3531             :             CALL helium_bead_solute_e_f(pint_env, helium, &
    3532        1027 :                                         endatom, ibead, pos(:, endatom, ibead), energy=energy)
    3533        1818 :             partaction = partaction + energy
    3534             :          END DO
    3535             :       END IF
    3536             : 
    3537        4090 :       partaction = partaction*helium%tau
    3538             : 
    3539        4090 :    END FUNCTION worm_path_inter_action_tail
    3540             : 
    3541             : END MODULE helium_worm

Generated by: LCOV version 1.15