LCOV - code coverage report
Current view: top level - src/motion - helium_worm.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 86.3 % 1571 1356
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 19 19

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

Generated by: LCOV version 2.0-1