LCOV - code coverage report
Current view: top level - src/motion - helium_common.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 623 966 64.5 %
Date: 2024-04-20 06:29:22 Functions: 23 40 57.5 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief  Independent helium subroutines shared with other modules
      10             : !> \author Lukasz Walewski
      11             : !> \date   2009-07-14
      12             : !> \note   Avoiding circular deps: do not USE any other helium_* modules here.
      13             : ! **************************************************************************************************
      14             : MODULE helium_common
      15             : 
      16             :    USE helium_types,                    ONLY: helium_solvent_type,&
      17             :                                               int_arr_ptr,&
      18             :                                               rho_atom_number,&
      19             :                                               rho_moment_of_inertia,&
      20             :                                               rho_num,&
      21             :                                               rho_projected_area,&
      22             :                                               rho_winding_cycle,&
      23             :                                               rho_winding_number
      24             :    USE input_constants,                 ONLY: denominator_inertia,&
      25             :                                               denominator_natoms,&
      26             :                                               denominator_rperp2,&
      27             :                                               helium_cell_shape_octahedron
      28             :    USE kinds,                           ONLY: default_string_length,&
      29             :                                               dp
      30             :    USE mathconstants,                   ONLY: pi
      31             :    USE memory_utilities,                ONLY: reallocate
      32             :    USE physcon,                         ONLY: angstrom,&
      33             :                                               bohr
      34             :    USE pint_public,                     ONLY: pint_com_pos
      35             :    USE pint_types,                      ONLY: pint_env_type
      36             :    USE splines_methods,                 ONLY: spline_value
      37             :    USE splines_types,                   ONLY: spline_data_type
      38             : #include "../base/base_uses.f90"
      39             : 
      40             :    IMPLICIT NONE
      41             : 
      42             :    PRIVATE
      43             : 
      44             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_common'
      46             : 
      47             :    PUBLIC :: helium_pbc
      48             :    PUBLIC :: helium_boxmean_3d
      49             :    PUBLIC :: helium_calc_rdf
      50             :    PUBLIC :: helium_calc_rho
      51             :    PUBLIC :: helium_calc_plength
      52             :    PUBLIC :: helium_rotate
      53             :    PUBLIC :: helium_eval_expansion
      54             :    PUBLIC :: helium_eval_chain
      55             :    PUBLIC :: helium_update_transition_matrix
      56             :    PUBLIC :: helium_update_coord_system
      57             :    PUBLIC :: helium_spline
      58             :    PUBLIC :: helium_cycle_number
      59             :    PUBLIC :: helium_path_length
      60             :    PUBLIC :: helium_is_winding
      61             :    PUBLIC :: helium_cycle_of
      62             :    PUBLIC :: helium_total_winding_number
      63             :    PUBLIC :: helium_total_projected_area
      64             :    PUBLIC :: helium_total_moment_of_inertia
      65             :    PUBLIC :: helium_com
      66             :    PUBLIC :: helium_set_rdf_coord_system
      67             : 
      68             : CONTAINS
      69             : 
      70             : ! ***************************************************************************
      71             : !> \brief  General PBC routine for helium.
      72             : !> \param helium ...
      73             : !> \param r ...
      74             : !> \param enforce ...
      75             : !> \date   2019-12-09
      76             : !> \author Harald Forbert
      77             : !> \note  Apply periodic boundary conditions as needed
      78             : !> \note  Combines several older routines into a single simpler/faster routine
      79             : ! **************************************************************************************************
      80   106955608 :    SUBROUTINE helium_pbc(helium, r, enforce)
      81             : 
      82             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
      83             :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT)         :: r
      84             :       LOGICAL, OPTIONAL                                  :: enforce
      85             : 
      86             :       REAL(KIND=dp)                                      :: cell_size, cell_size_inv, corr, rx, ry, &
      87             :                                                             rz, sx, sy, sz
      88             : 
      89   106955608 :       IF (.NOT. (helium%periodic .OR. PRESENT(enforce))) RETURN
      90             : 
      91   101673513 :       cell_size = helium%cell_size
      92   101673513 :       cell_size_inv = helium%cell_size_inv
      93             : 
      94   101673513 :       rx = r(1)*cell_size_inv
      95   101673513 :       IF (rx > 0.5_dp) THEN
      96     4098339 :          rx = rx - REAL(INT(rx + 0.5_dp), dp)
      97    97575174 :       ELSEIF (rx < -0.5_dp) THEN
      98     4159268 :          rx = rx - REAL(INT(rx - 0.5_dp), dp)
      99             :       END IF
     100             : 
     101   101673513 :       ry = r(2)*cell_size_inv
     102   101673513 :       IF (ry > 0.5_dp) THEN
     103     3960982 :          ry = ry - REAL(INT(ry + 0.5_dp), dp)
     104    97712531 :       ELSEIF (ry < -0.5_dp) THEN
     105     8203072 :          ry = ry - REAL(INT(ry - 0.5_dp), dp)
     106             :       END IF
     107             : 
     108   101673513 :       rz = r(3)*cell_size_inv
     109   101673513 :       IF (rz > 0.5_dp) THEN
     110     4236780 :          rz = rz - REAL(INT(rz + 0.5_dp), dp)
     111    97436733 :       ELSEIF (rz < -0.5_dp) THEN
     112    10562463 :          rz = rz - REAL(INT(rz - 0.5_dp), dp)
     113             :       END IF
     114             : 
     115   101673513 :       IF (helium%cell_shape == helium_cell_shape_octahedron) THEN
     116   101673513 :          corr = 0.0_dp
     117   101673513 :          IF (rx > 0.0_dp) THEN
     118    50015863 :             corr = corr + rx
     119    50015863 :             sx = 0.5_dp
     120             :          ELSE
     121    51657650 :             corr = corr - rx
     122    51657650 :             sx = -0.5_dp
     123             :          END IF
     124   101673513 :          IF (ry > 0.0_dp) THEN
     125    48620754 :             corr = corr + ry
     126    48620754 :             sy = 0.5_dp
     127             :          ELSE
     128    53052759 :             corr = corr - ry
     129    53052759 :             sy = -0.5_dp
     130             :          END IF
     131   101673513 :          IF (rz > 0.0_dp) THEN
     132    50109946 :             corr = corr + rz
     133    50109946 :             sz = 0.5_dp
     134             :          ELSE
     135    51563567 :             corr = corr - rz
     136    51563567 :             sz = -0.5_dp
     137             :          END IF
     138   101673513 :          IF (corr > 0.75_dp) THEN
     139    36034635 :             rx = rx - sx
     140    36034635 :             ry = ry - sy
     141    36034635 :             rz = rz - sz
     142             :          END IF
     143             :       END IF
     144             : 
     145   101673513 :       r(1) = rx*cell_size
     146   101673513 :       r(2) = ry*cell_size
     147   101673513 :       r(3) = rz*cell_size
     148             : 
     149             :    END SUBROUTINE helium_pbc
     150             : 
     151             : ! ***************************************************************************
     152             : !> \brief  find distance squared of nearest image
     153             : !> \param helium ...
     154             : !> \param r ...
     155             : !> \return ...
     156             : !> \date   2019-12-09
     157             : !> \author Harald Forbert
     158             : !> \note   Given a distance vector r, return the distance squared
     159             : !>         of the nearest image. Cell information is taken from the
     160             : !>         helium environment argument.
     161             : ! **************************************************************************************************
     162             : 
     163   367358409 :    PURE FUNCTION helium_pbc_r2(helium, r)
     164             : 
     165             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     166             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: r
     167             :       REAL(KIND=dp)                                      :: helium_pbc_r2
     168             : 
     169             :       REAL(KIND=dp)                                      :: cell_size_inv, corr, rx, ry, rz
     170             : 
     171   367358409 :       IF (helium%periodic) THEN
     172   348193214 :          cell_size_inv = helium%cell_size_inv
     173   348193214 :          rx = ABS(r(1)*cell_size_inv)
     174   348193214 :          rx = ABS(rx - REAL(INT(rx + 0.5_dp), dp))
     175   348193214 :          ry = ABS(r(2)*cell_size_inv)
     176   348193214 :          ry = ABS(ry - REAL(INT(ry + 0.5_dp), dp))
     177   348193214 :          rz = ABS(r(3)*cell_size_inv)
     178   348193214 :          rz = ABS(rz - REAL(INT(rz + 0.5_dp), dp))
     179   348193214 :          IF (helium%cell_shape == helium_cell_shape_octahedron) THEN
     180   348193214 :             corr = rx + ry + rz - 0.75_dp
     181   348193214 :             IF (corr < 0.0_dp) corr = 0.0_dp
     182   348193214 :             helium_pbc_r2 = helium%cell_size**2*(rx*rx + ry*ry + rz*rz - corr)
     183             :          ELSE
     184           0 :             helium_pbc_r2 = helium%cell_size**2*(rx*rx + ry*ry + rz*rz)
     185             :          END IF
     186             :       ELSE
     187    19165195 :          helium_pbc_r2 = (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     188             :       END IF
     189   367358409 :    END FUNCTION helium_pbc_r2
     190             : 
     191             : ! ***************************************************************************
     192             : !> \brief  Calculate the point equidistant from two other points a and b
     193             : !>         within the helium box - 3D version
     194             : !> \param    helium   helium environment for which
     195             : !> \param a vectors for which to find the mean within the He box
     196             : !> \param b vectors for which to find the mean within the He box
     197             : !> \param c ...
     198             : !> \par History
     199             : !>      2009-10-02 renamed, originally was helium_boxmean [lwalewski]
     200             : !>      2009-10-02 redesigned so it is now called as a subroutine [lwalewski]
     201             : !>      2009-10-02 redesigned so it now gets/returns a 3D vectors [lwalewski]
     202             : !> \author hforbert
     203             : ! **************************************************************************************************
     204     2959734 :    SUBROUTINE helium_boxmean_3d(helium, a, b, c)
     205             : 
     206             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     207             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a, b
     208             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: c
     209             : 
     210    11838936 :       c(:) = b(:) - a(:)
     211     2959734 :       CALL helium_pbc(helium, c)
     212    11838936 :       c(:) = a(:) + 0.5_dp*c(:)
     213     2959734 :       CALL helium_pbc(helium, c)
     214     2959734 :    END SUBROUTINE helium_boxmean_3d
     215             : 
     216             : ! ***************************************************************************
     217             : !> \brief  Given the permutation state assign cycle lengths to all He atoms.
     218             : !> \param helium ...
     219             : !> \date   2011-07-06
     220             : !> \author Lukasz Walewski
     221             : !> \note   The helium%atom_plength array is filled with cycle lengths,
     222             : !>         each atom gets the length of the permutation cycle it belongs to.
     223             : ! **************************************************************************************************
     224          12 :    SUBROUTINE helium_calc_atom_cycle_length(helium)
     225             : 
     226             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     227             : 
     228             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_atom_cycle_length'
     229             : 
     230             :       CHARACTER(len=default_string_length)               :: err_str
     231             :       INTEGER                                            :: clen, curr_idx, handle, ia, start_idx
     232          12 :       INTEGER, DIMENSION(:), POINTER                     :: atoms_in_cycle
     233             :       LOGICAL                                            :: atoms_left, path_end_reached
     234          12 :       LOGICAL, DIMENSION(:), POINTER                     :: atom_was_used
     235             : 
     236          12 :       CALL timeset(routineN, handle)
     237             : 
     238          12 :       NULLIFY (atoms_in_cycle)
     239          12 :       atoms_in_cycle => helium%itmp_atoms_1d
     240         396 :       atoms_in_cycle(:) = 0
     241             : 
     242          12 :       NULLIFY (atom_was_used)
     243          12 :       atom_was_used => helium%ltmp_atoms_1d
     244         396 :       atom_was_used(:) = .FALSE.
     245             : 
     246         396 :       helium%atom_plength(:) = 0
     247             : 
     248             :       start_idx = 1
     249             :       DO
     250         384 :          clen = 0
     251         384 :          path_end_reached = .FALSE.
     252         384 :          curr_idx = start_idx
     253         384 :          DO ia = 1, helium%atoms
     254         384 :             clen = clen + 1
     255         384 :             atoms_in_cycle(clen) = curr_idx
     256         384 :             atom_was_used(curr_idx) = .TRUE.
     257             : 
     258             :             ! follow to the next atom in the cycle
     259         384 :             curr_idx = helium%permutation(curr_idx)
     260         384 :             IF (curr_idx .EQ. start_idx) THEN
     261             :                path_end_reached = .TRUE.
     262             :                EXIT
     263             :             END IF
     264             :          END DO
     265         384 :          err_str = "Permutation path is not cyclic."
     266         384 :          IF (.NOT. path_end_reached) THEN
     267           0 :             CPABORT(err_str)
     268             :          END IF
     269             : 
     270             :          ! assign the cycle length to the collected atoms
     271         768 :          DO ia = 1, clen
     272         768 :             helium%atom_plength(atoms_in_cycle(ia)) = clen
     273             :          END DO
     274             : 
     275             :          ! go to the next "unused" atom
     276         384 :          atoms_left = .FALSE.
     277        6720 :          DO ia = 1, helium%atoms
     278        6720 :             IF (.NOT. atom_was_used(ia)) THEN
     279             :                start_idx = ia
     280             :                atoms_left = .TRUE.
     281             :                EXIT
     282             :             END IF
     283             :          END DO
     284             : 
     285         384 :          IF (.NOT. atoms_left) EXIT
     286             :       END DO
     287             : 
     288         396 :       atoms_in_cycle(:) = 0
     289         396 :       NULLIFY (atoms_in_cycle)
     290         396 :       atom_was_used(:) = .FALSE.
     291          12 :       NULLIFY (atom_was_used)
     292             : 
     293          12 :       CALL timestop(handle)
     294             : 
     295          12 :    END SUBROUTINE helium_calc_atom_cycle_length
     296             : 
     297             : ! ***************************************************************************
     298             : !> \brief  Decompose a permutation into cycles
     299             : !> \param permutation ...
     300             : !> \return ...
     301             : !> \date   2013-12-11
     302             : !> \author Lukasz Walewski
     303             : !> \note   Given a permutation return a pointer to an array of pointers,
     304             : !>         with each element pointing to a cycle of this permutation.
     305             : !> \note   This function allocates memory and returns a pointer to it,
     306             : !>         do not forget to deallocate once finished with the results.
     307             : ! **************************************************************************************************
     308           0 :    FUNCTION helium_calc_cycles(permutation) RESULT(cycles)
     309             : 
     310             :       INTEGER, DIMENSION(:), POINTER                     :: permutation
     311             :       TYPE(int_arr_ptr), DIMENSION(:), POINTER           :: cycles
     312             : 
     313             :       INTEGER                                            :: curat, ncycl, nsize, nused
     314           0 :       INTEGER, DIMENSION(:), POINTER                     :: cur_cycle, used_indices
     315           0 :       TYPE(int_arr_ptr), DIMENSION(:), POINTER           :: my_cycles
     316             : 
     317           0 :       nsize = SIZE(permutation)
     318             : 
     319             :       ! most pesimistic scenario: no cycles longer than 1
     320           0 :       ALLOCATE (my_cycles(nsize))
     321             : 
     322             :       ! go over the permutation and identify cycles
     323           0 :       curat = 1
     324           0 :       nused = 0
     325           0 :       ncycl = 0
     326           0 :       DO WHILE (curat .LE. nsize)
     327             : 
     328             :          ! get the permutation cycle the current atom belongs to
     329           0 :          cur_cycle => helium_cycle_of(curat, permutation)
     330             : 
     331             :          ! include the current cycle in the pool of "used" indices
     332           0 :          nused = nused + SIZE(cur_cycle)
     333           0 :          CALL reallocate(used_indices, 1, nused)
     334           0 :          used_indices(nused - SIZE(cur_cycle) + 1:nused) = cur_cycle(1:SIZE(cur_cycle))
     335             : 
     336             :          ! store the pointer to the current cycle
     337           0 :          ncycl = ncycl + 1
     338           0 :          my_cycles(ncycl)%iap => cur_cycle
     339             : 
     340             :          ! free the local pointer
     341           0 :          NULLIFY (cur_cycle)
     342             : 
     343             :          ! try to increment the current atom index
     344           0 :          DO WHILE (ANY(used_indices .EQ. curat))
     345           0 :             curat = curat + 1
     346             :          END DO
     347             : 
     348             :       END DO
     349             : 
     350           0 :       DEALLOCATE (used_indices)
     351             :       NULLIFY (used_indices)
     352             : 
     353             :       ! assign the result
     354           0 :       ALLOCATE (cycles(ncycl))
     355           0 :       cycles(1:ncycl) = my_cycles(1:ncycl)
     356             : 
     357           0 :       DEALLOCATE (my_cycles)
     358           0 :       NULLIFY (my_cycles)
     359             : 
     360           0 :    END FUNCTION helium_calc_cycles
     361             : 
     362             : ! ***************************************************************************
     363             : !> \brief  Calculate helium density distribution functions and store them
     364             : !>         in helium%rho_inst
     365             : !> \param helium ...
     366             : !> \date   2011-06-14
     367             : !> \author Lukasz Walewski
     368             : ! **************************************************************************************************
     369          12 :    SUBROUTINE helium_calc_rho(helium)
     370             : 
     371             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     372             : 
     373             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_calc_rho'
     374             : 
     375             :       CHARACTER(len=default_string_length)               :: msgstr
     376             :       INTEGER                                            :: aa, bx, by, bz, handle, ia, ib, ic, id, &
     377             :                                                             ii, ir, n_out_of_range, nbin
     378             :       INTEGER, DIMENSION(3)                              :: nw
     379          12 :       INTEGER, DIMENSION(:), POINTER                     :: cycl_len
     380             :       LOGICAL                                            :: ltmp1, ltmp2, ltmp3
     381             :       REAL(KIND=dp)                                      :: invd3r, invdr, invp, rtmp
     382             :       REAL(KIND=dp), DIMENSION(3)                        :: maxr_half, r, vlink, vtotal, wn
     383          12 :       TYPE(int_arr_ptr), DIMENSION(:), POINTER           :: cycles
     384             : 
     385          12 :       CALL timeset(routineN, handle)
     386             : 
     387          48 :       maxr_half(:) = helium%rho_maxr/2.0_dp
     388          12 :       invdr = 1.0_dp/helium%rho_delr
     389          12 :       invd3r = invdr*invdr*invdr
     390          12 :       invp = 1.0_dp/REAL(helium%beads, dp)
     391          12 :       nbin = helium%rho_nbin
     392             :       ! assign the cycle lengths to all the atoms
     393          12 :       CALL helium_calc_atom_cycle_length(helium)
     394             :       NULLIFY (cycl_len)
     395          12 :       cycl_len => helium%atom_plength
     396          72 :       DO ir = 1, rho_num ! loop over densities ---
     397             : 
     398          60 :          IF (.NOT. helium%rho_property(ir)%is_calculated) THEN
     399             :             ! skip densities that are not requested by the user
     400             :             CYCLE
     401             :          END IF
     402             : 
     403          12 :          SELECT CASE (ir)
     404             : 
     405             :          CASE (rho_atom_number)
     406          12 :             ii = helium%rho_property(ir)%component_index(1)
     407        9912 :             helium%rho_incr(ii, :, :) = invp
     408             : 
     409             :          CASE (rho_projected_area)
     410           0 :             vtotal(:) = helium_total_projected_area(helium)
     411           0 :             DO ia = 1, helium%atoms
     412           0 :                DO ib = 1, helium%beads
     413           0 :                   vlink(:) = helium_link_projected_area(helium, ia, ib)
     414           0 :                   DO ic = 1, 3
     415           0 :                      ii = helium%rho_property(ir)%component_index(ic)
     416           0 :                      helium%rho_incr(ii, ia, ib) = vtotal(ic)*vlink(ic)*angstrom*angstrom*angstrom*angstrom
     417             :                   END DO
     418             :                END DO
     419             :             END DO
     420             : 
     421             : !        CASE (rho_winding_number)
     422             : !          vtotal(:) = helium_total_winding_number(helium)
     423             : !          DO ia = 1, helium%atoms
     424             : !            DO ib = 1, helium%beads
     425             : !              vlink(:) = helium_link_winding_number(helium,ia,ib)
     426             : !              DO ic = 1, 3
     427             : !                ii = helium%rho_property(ir)%component_index(ic)
     428             : !                helium%rho_incr(ii,ia,ib) = vtotal(ic)*vlink(ic)*angstrom*angstrom
     429             : !              END DO
     430             : !            END DO
     431             : !          END DO
     432             : 
     433             :          CASE (rho_winding_number)
     434           0 :             vtotal(:) = helium_total_winding_number(helium)
     435           0 :             DO id = 1, 3
     436           0 :                ii = helium%rho_property(ir)%component_index(id)
     437           0 :                helium%rho_incr(ii, :, :) = 0.0_dp
     438             :             END DO
     439           0 :             NULLIFY (cycles)
     440           0 :             cycles => helium_calc_cycles(helium%permutation)
     441           0 :             DO ic = 1, SIZE(cycles)
     442           0 :                wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
     443           0 :                DO ia = 1, SIZE(cycles(ic)%iap)
     444           0 :                   aa = cycles(ic)%iap(ia)
     445           0 :                   DO ib = 1, helium%beads
     446           0 :                      vlink(:) = helium_link_winding_number(helium, aa, ib)
     447           0 :                      DO id = 1, 3
     448           0 :                         IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
     449           0 :                            ii = helium%rho_property(ir)%component_index(id)
     450           0 :                            helium%rho_incr(ii, aa, ib) = vtotal(id)*vlink(id)*angstrom*angstrom
     451             :                         END IF
     452             :                      END DO
     453             :                   END DO
     454             :                END DO
     455             :             END DO
     456           0 :             DEALLOCATE (cycles)
     457             : 
     458             :          CASE (rho_winding_cycle)
     459           0 :             vtotal(:) = helium_total_winding_number(helium)
     460           0 :             DO id = 1, 3
     461           0 :                ii = helium%rho_property(ir)%component_index(id)
     462           0 :                helium%rho_incr(ii, :, :) = 0.0_dp
     463             :             END DO
     464           0 :             NULLIFY (cycles)
     465           0 :             cycles => helium_calc_cycles(helium%permutation)
     466             :             ! compute number of atoms in all winding cycles
     467           0 :             nw(:) = 0
     468           0 :             DO ic = 1, SIZE(cycles)
     469           0 :                wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
     470           0 :                DO id = 1, 3
     471           0 :                   IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
     472           0 :                      nw(id) = nw(id) + SIZE(cycles(ic)%iap)
     473             :                   END IF
     474             :                END DO
     475             :             END DO
     476             :             ! assign contributions to all beads of winding cycles only
     477           0 :             DO ic = 1, SIZE(cycles)
     478           0 :                wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
     479           0 :                DO id = 1, 3
     480           0 :                   IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
     481           0 :                      DO ia = 1, SIZE(cycles(ic)%iap)
     482           0 :                         aa = cycles(ic)%iap(ia)
     483           0 :                         DO ib = 1, helium%beads
     484           0 :                            IF (nw(id) .GT. 0) THEN ! this test should always get passed
     485           0 :                               ii = helium%rho_property(ir)%component_index(id)
     486           0 :                               rtmp = invp/REAL(nw(id), dp)
     487           0 :                               helium%rho_incr(ii, aa, ib) = rtmp*vtotal(id)*vtotal(id)*angstrom*angstrom
     488             :                            END IF
     489             :                         END DO
     490             :                      END DO
     491             :                   END IF
     492             :                END DO
     493             :             END DO
     494           0 :             DEALLOCATE (cycles)
     495             : 
     496             :          CASE (rho_moment_of_inertia)
     497           0 :             vtotal(:) = helium_total_moment_of_inertia(helium)
     498          12 :             DO ia = 1, helium%atoms
     499           0 :                DO ib = 1, helium%beads
     500           0 :                   vlink(:) = helium_link_moment_of_inertia(helium, ia, ib)
     501           0 :                   DO ic = 1, 3
     502           0 :                      ii = helium%rho_property(ir)%component_index(ic)
     503           0 :                      helium%rho_incr(ii, ia, ib) = vlink(ic)*angstrom*angstrom
     504             :                   END DO
     505             :                END DO
     506             :             END DO
     507             : 
     508             :          CASE DEFAULT
     509             :             ! do nothing
     510             : 
     511             :          END SELECT
     512             : 
     513             :       END DO ! loop over densities ---
     514             : 
     515          12 :       n_out_of_range = 0
     516       25332 :       helium%rho_inst(:, :, :, :) = 0.0_dp
     517         396 :       DO ia = 1, helium%atoms
     518             :          ! bin the bead positions of the current atom using the increments set above
     519        9996 :          DO ib = 1, helium%beads
     520             :             ! map the current bead position to the corresponding voxel
     521       38400 :             r(:) = helium%pos(:, ia, ib) - helium%center(:)
     522             :             ! enforce PBC even if this is a non-periodic calc to avoid density leakage
     523        9600 :             CALL helium_pbc(helium, r, enforce=.TRUE.)
     524             :             ! set up bin indices (translate by L/2 to avoid non-positive array indices)
     525        9600 :             bx = INT((r(1) + maxr_half(1))*invdr) + 1
     526        9600 :             by = INT((r(2) + maxr_half(2))*invdr) + 1
     527        9600 :             bz = INT((r(3) + maxr_half(3))*invdr) + 1
     528             :             ! check that the resulting bin numbers are within array bounds
     529        9600 :             ltmp1 = (0 .LT. bx) .AND. (bx .LE. nbin)
     530        9600 :             ltmp2 = (0 .LT. by) .AND. (by .LE. nbin)
     531        9600 :             ltmp3 = (0 .LT. bz) .AND. (bz .LE. nbin)
     532        9984 :             IF (ltmp1 .AND. ltmp2 .AND. ltmp3) THEN
     533             :                ! increment all the estimators (those that the current atom does not
     534             :                ! contribute to have increment incr(ic)==0)
     535       19200 :                DO ic = 1, helium%rho_num_act
     536       19200 :                   helium%rho_inst(ic, bx, by, bz) = helium%rho_inst(ic, bx, by, bz) + helium%rho_incr(ic, ia, ib)
     537             :                END DO
     538             :             ELSE
     539           0 :                n_out_of_range = n_out_of_range + 1
     540             :             END IF
     541             :          END DO
     542             :       END DO
     543             :       ! scale by volume element
     544       25332 :       helium%rho_inst(:, :, :, :) = helium%rho_inst(:, :, :, :)*invd3r
     545             : 
     546             :       ! stop if any bead fell out of the range
     547             :       ! since enforced PBC should have taken care of such leaks
     548          12 :       WRITE (msgstr, *) n_out_of_range
     549          12 :       msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr))
     550          12 :       IF (n_out_of_range .GT. 0) THEN
     551           0 :          CPABORT(msgstr)
     552             :       END IF
     553             : 
     554          12 :       CALL timestop(handle)
     555             : 
     556          12 :    END SUBROUTINE helium_calc_rho
     557             : 
     558             : #if 0
     559             : ! ***************************************************************************
     560             : !> \brief  Normalize superfluid densities according to the input keyword
     561             : !>         HELIUM%SUPERFLUID_ESTIMATOR%DENOMINATOR
     562             : !> \param helium ...
     563             : !> \param rho ...
     564             : !> \date   2014-06-24
     565             : !> \author Lukasz Walewski
     566             : ! **************************************************************************************************
     567             :    SUBROUTINE helium_norm_rho(helium, rho)
     568             : 
     569             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     570             :       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: rho
     571             : 
     572             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_norm_rho', &
     573             :          routineP = moduleN//':'//routineN
     574             : 
     575             :       INTEGER                                            :: ix, iy, iz, ndim
     576             :       REAL(KIND=dp)                                      :: invatoms, rx, ry, rz
     577             :       REAL(KIND=dp), DIMENSION(3)                        :: invmoit, invrperp, ro
     578             : 
     579             :       SELECT CASE (helium%supest_denominator)
     580             : 
     581             :       CASE (denominator_natoms)
     582             :          invatoms = 1.0_dp/REAL(helium%atoms, dp)
     583             :          rho(2, :, :, :) = rho(2, :, :, :)*invatoms
     584             :          rho(3, :, :, :) = rho(3, :, :, :)*invatoms
     585             :          rho(4, :, :, :) = rho(4, :, :, :)*invatoms
     586             : 
     587             :       CASE (denominator_inertia)
     588             :          invmoit(:) = REAL(helium%atoms, dp)/helium%mominer%ravr(:)
     589             :          rho(2, :, :, :) = rho(2, :, :, :)*invmoit(1)
     590             :          rho(3, :, :, :) = rho(3, :, :, :)*invmoit(2)
     591             :          rho(4, :, :, :) = rho(4, :, :, :)*invmoit(3)
     592             : 
     593             :       CASE (denominator_rperp2)
     594             :          ndim = helium%rho_nbin
     595             :          ro(:) = helium%center(:) - 0.5_dp*(helium%rho_maxr - helium%rho_delr)
     596             :          DO ix = 1, ndim
     597             :             DO iy = 1, ndim
     598             :                DO iz = 1, ndim
     599             :                   rx = ro(1) + REAL(ix - 1, dp)*helium%rho_delr
     600             :                   ry = ro(2) + REAL(iy - 1, dp)*helium%rho_delr
     601             :                   rz = ro(3) + REAL(iz - 1, dp)*helium%rho_delr
     602             :                   invrperp(1) = 1.0_dp/(ry*ry + rz*rz)
     603             :                   invrperp(2) = 1.0_dp/(rz*rz + rx*rx)
     604             :                   invrperp(3) = 1.0_dp/(rx*rx + ry*ry)
     605             :                   rho(2, ix, iy, iz) = rho(2, ix, iy, iz)*invrperp(1)
     606             :                   rho(3, ix, iy, iz) = rho(3, ix, iy, iz)*invrperp(2)
     607             :                   rho(4, ix, iy, iz) = rho(4, ix, iy, iz)*invrperp(3)
     608             :                END DO
     609             :             END DO
     610             :          END DO
     611             : 
     612             :       CASE DEFAULT
     613             :          ! do nothing
     614             : 
     615             :       END SELECT
     616             : 
     617             :    END SUBROUTINE helium_norm_rho
     618             : #endif
     619             : 
     620             : ! ***************************************************************************
     621             : !> \brief  Calculate helium radial distribution functions wrt positions given
     622             : !>         by <centers> using the atomic weights given by <weights>.
     623             : !> \param helium ...
     624             : !> \date   2009-07-22
     625             : !> \par    History
     626             : !>         2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
     627             : !> \author Lukasz Walewski
     628             : ! **************************************************************************************************
     629          12 :    SUBROUTINE helium_calc_rdf(helium)
     630             : 
     631             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     632             : 
     633             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_calc_rdf'
     634             : 
     635             :       CHARACTER(len=default_string_length)               :: msgstr
     636             :       INTEGER                                            :: bin, handle, ia, ib, ic, ind_hehe, &
     637             :                                                             n_out_of_range, nbin
     638             :       REAL(KIND=dp)                                      :: invdr, nideal, ri, rlower, rupper
     639          12 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: pref
     640             :       REAL(KIND=dp), DIMENSION(3)                        :: r, r0
     641             : 
     642          12 :       CALL timeset(routineN, handle)
     643             : 
     644          12 :       invdr = 1.0_dp/helium%rdf_delr
     645          12 :       nbin = helium%rdf_nbin
     646          12 :       n_out_of_range = 0
     647        6012 :       helium%rdf_inst(:, :) = 0.0_dp
     648             : 
     649          12 :       ind_hehe = 0
     650             :       ! Calculate He-He RDF
     651          12 :       IF (helium%rdf_he_he) THEN
     652          12 :          ind_hehe = 1
     653         312 :          DO ib = 1, helium%beads
     654        9912 :             DO ia = 1, helium%atoms
     655       38400 :                r0(:) = helium%pos(:, ia, ib)
     656             : 
     657      317100 :                DO ic = 1, helium%atoms
     658      307200 :                   IF (ia == ic) CYCLE
     659     1190400 :                   r(:) = helium%pos(:, ic, ib) - r0(:)
     660      297600 :                   CALL helium_pbc(helium, r)
     661      297600 :                   ri = SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     662      297600 :                   bin = INT(ri*invdr) + 1
     663      307200 :                   IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
     664             :                      ! increment the RDF value for He atoms inside the r_6 sphere
     665      194384 :                      helium%rdf_inst(ind_hehe, bin) = helium%rdf_inst(ind_hehe, bin) + 1.0_dp
     666             :                   ELSE
     667      103216 :                      n_out_of_range = n_out_of_range + 1
     668             :                   END IF
     669             :                END DO
     670             :             END DO
     671             :          END DO
     672             :       END IF
     673             : 
     674             :       ! Calculate Solute-He RDF
     675          12 :       IF (helium%solute_present .AND. helium%rdf_sol_he) THEN
     676           0 :          DO ib = 1, helium%beads
     677           0 :             DO ia = 1, helium%solute_atoms
     678           0 :                r0(:) = helium%rdf_centers(ib, 3*(ia - 1) + 1:3*(ia - 1) + 3)
     679             : 
     680           0 :                DO ic = 1, helium%atoms
     681           0 :                   r(:) = helium%pos(:, ic, ib) - r0(:)
     682           0 :                   CALL helium_pbc(helium, r)
     683           0 :                   ri = SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     684           0 :                   bin = INT(ri*invdr) + 1
     685           0 :                   IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
     686             :                      ! increment the RDF value for He atoms inside the r_6 sphere
     687           0 :                      helium%rdf_inst(ind_hehe + ia, bin) = helium%rdf_inst(ind_hehe + ia, bin) + 1.0_dp
     688             :                   ELSE
     689           0 :                      n_out_of_range = n_out_of_range + 1
     690             :                   END IF
     691             :                END DO
     692             :             END DO
     693             :          END DO
     694             : 
     695             :       END IF
     696             : 
     697             :       ! for periodic case we intentionally truncate RDFs to spherical volume
     698             :       ! so we skip atoms in the corners
     699          12 :       IF (.NOT. helium%periodic) THEN
     700           0 :          IF (n_out_of_range .GT. 0) THEN
     701           0 :             WRITE (msgstr, *) n_out_of_range
     702           0 :             msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr))
     703           0 :             CPWARN(msgstr)
     704             :          END IF
     705             :       END IF
     706             : 
     707             :       ! normalize the histogram to get g(r)
     708          36 :       ALLOCATE (pref(helium%rdf_num))
     709          24 :       pref(:) = 0.0_dp
     710          12 :       IF (helium%periodic) THEN
     711             :          ! Use helium density for normalization
     712          24 :          pref(:) = helium%density*helium%beads*helium%atoms
     713             :          ! Correct for He-He-RDF
     714          12 :          IF (helium%rdf_he_he) THEN
     715          12 :             pref(1) = pref(1)/helium%atoms*MAX(helium%atoms - 1, 1)
     716             :          END IF
     717             :       ELSE
     718             :          ! Non-periodic case has density of 0, use integral for normalzation
     719             :          ! This leads to a unit of 1/volume of the RDF
     720           0 :          pref(:) = 0.5_dp*helium%rdf_inst(:, 1)
     721           0 :          DO bin = 2, helium%rdf_nbin - 1
     722           0 :             pref(:) = pref(:) + helium%rdf_inst(:, bin)
     723             :          END DO
     724           0 :          pref(:) = pref(:) + 0.5_dp*helium%rdf_inst(:, helium%rdf_nbin)
     725             : 
     726             :          !set integral of histogram to number of atoms:
     727           0 :          pref(:) = pref(:)/helium%atoms
     728             :          ! Correct for He-He-RDF
     729           0 :          IF (helium%rdf_he_he) THEN
     730           0 :             pref(1) = pref(1)*helium%atoms/MAX(helium%atoms - 1, 1)
     731             :          END IF
     732             :       END IF
     733             :       ! Volume integral first:
     734        3012 :       DO bin = 1, helium%rdf_nbin
     735        3000 :          rlower = REAL(bin - 1, dp)*helium%rdf_delr
     736        3000 :          rupper = rlower + helium%rdf_delr
     737        3000 :          nideal = (rupper**3 - rlower**3)*4.0_dp*pi/3.0_dp
     738        6012 :          helium%rdf_inst(:, bin) = helium%rdf_inst(:, bin)/nideal
     739             :       END DO
     740             :       ! No normalization for density
     741          24 :       pref(:) = 1.0_dp/pref(:)
     742          24 :       DO ia = 1, helium%rdf_num
     743        3024 :          helium%rdf_inst(ia, :) = helium%rdf_inst(ia, :)*pref(ia)
     744             :       END DO
     745             : 
     746          12 :       DEALLOCATE (pref)
     747             : 
     748          12 :       CALL timestop(handle)
     749             : 
     750          12 :    END SUBROUTINE helium_calc_rdf
     751             : 
     752             : ! ***************************************************************************
     753             : !> \brief  Calculate probability distribution of the permutation lengths
     754             : !> \param helium ...
     755             : !> \date   2010-06-07
     756             : !> \author Lukasz Walewski
     757             : !> \note   Valid permutation path length is an integer (1, NATOMS), number
     758             : !>         of paths of a given length is calculated here and average over
     759             : !>         inner loop iterations and helium environments is done in
     760             : !>         helium_sample.
     761             : ! **************************************************************************************************
     762        7047 :    PURE SUBROUTINE helium_calc_plength(helium)
     763             : 
     764             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     765             : 
     766             :       INTEGER                                            :: i, j, k
     767             : 
     768      221151 :       helium%plength_inst(:) = 0.0_dp
     769      221151 :       DO i = 1, helium%atoms
     770      214104 :          j = helium%permutation(i)
     771      214104 :          k = 1
     772           0 :          DO
     773      214104 :             IF (j == i) EXIT
     774           0 :             k = k + 1
     775           0 :             j = helium%permutation(j)
     776             :          END DO
     777      221151 :          helium%plength_inst(k) = helium%plength_inst(k) + 1
     778             :       END DO
     779      221151 :       helium%plength_inst(:) = helium%plength_inst(:)/helium%atoms
     780             : 
     781        7047 :    END SUBROUTINE helium_calc_plength
     782             : 
     783             : ! ***************************************************************************
     784             : !> \brief  Rotate helium particles in imaginary time by nslices
     785             : !> \param helium ...
     786             : !> \param nslices ...
     787             : !> \author hforbert
     788             : !> \note   Positions of helium beads in helium%pos array are reorganized such
     789             : !>         that the indices are cyclically translated in a permutation-aware
     790             : !>         manner. helium%relrot is given a new value that represents the new
     791             : !>         'angle' of the beads. This is done modulo helium%beads, so relrot
     792             : !>         should be always within 0 (no rotation) and helium%beads-1 (almost
     793             : !>         full rotation). [lwalewski]
     794             : ! **************************************************************************************************
     795        7146 :    PURE SUBROUTINE helium_rotate(helium, nslices)
     796             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     797             :       INTEGER, INTENT(IN)                                :: nslices
     798             : 
     799             :       INTEGER                                            :: b, i, j, k, n
     800             : 
     801        7146 :       b = helium%beads
     802        7146 :       n = helium%atoms
     803        7146 :       i = MOD(nslices, b)
     804        7146 :       IF (i < 0) i = i + b
     805        7146 :       IF ((i >= b) .OR. (i < 1)) RETURN
     806        6750 :       helium%relrot = MOD(helium%relrot + i, b)
     807       75215 :       DO k = 1, i
     808     8470671 :          helium%work(:, :, k) = helium%pos(:, :, k)
     809             :       END DO
     810       75850 :       DO k = i + 1, b
     811     8514570 :          helium%pos(:, :, k - i) = helium%pos(:, :, k)
     812             :       END DO
     813       75215 :       DO k = 1, i
     814     2174079 :          DO j = 1, n
     815     8463921 :             helium%pos(:, j, b - i + k) = helium%work(:, helium%permutation(j), k)
     816             :          END DO
     817             :       END DO
     818             :    END SUBROUTINE helium_rotate
     819             : 
     820             : ! ***************************************************************************
     821             : !> \brief  Calculate the pair-product action or energy by evaluating the
     822             : !>         power series expansion according to Eq. 4.46 in Ceperley 1995.
     823             : !> \param helium ...
     824             : !> \param r ...
     825             : !> \param rp ...
     826             : !> \param work ...
     827             : !> \param action ...
     828             : !> \return ...
     829             : !> \author Harald Forbert
     830             : ! **************************************************************************************************
     831    73517397 :    FUNCTION helium_eval_expansion(helium, r, rp, work, action) RESULT(res)
     832             : 
     833             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     834             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: r, rp
     835             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:), &
     836             :          INTENT(INOUT)                                   :: work
     837             :       LOGICAL, INTENT(IN), OPTIONAL                      :: action
     838             :       REAL(KIND=dp)                                      :: res
     839             : 
     840             :       INTEGER                                            :: i, j, nsp, tline
     841             :       LOGICAL                                            :: nocut
     842             :       REAL(KIND=dp)                                      :: a, ar, arp, b, h26, q, s, v, z
     843             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     844    73517397 :          POINTER                                         :: offdiagsplines
     845             :       TYPE(spline_data_type), POINTER                    :: endpspline
     846             : 
     847    73517397 :       endpspline => helium%u0
     848    73517397 :       offdiagsplines => helium%uoffdiag
     849    73517397 :       nocut = .TRUE.
     850    73517397 :       IF (PRESENT(action)) THEN
     851           0 :          IF (.NOT. action) THEN
     852           0 :             endpspline => helium%e0
     853           0 :             offdiagsplines => helium%eoffdiag
     854           0 :             nocut = .FALSE.
     855             :          END IF
     856             :       END IF
     857             : 
     858    73517397 :       ar = SQRT(helium_pbc_r2(helium, r))
     859    73517397 :       arp = SQRT(helium_pbc_r2(helium, rp))
     860             : 
     861    73517397 :       IF (helium%periodic .AND. ((ar > 0.5_dp*helium%cell_size) &
     862             :                                  .OR. (arp > 0.5_dp*helium%cell_size))) THEN
     863     6754137 :          v = 0.0_dp
     864     6754137 :          IF (arp > 0.5_dp*helium%cell_size) THEN
     865     4649159 :             IF (nocut) v = v + helium_spline(endpspline, 0.5_dp*helium%cell_size)
     866             :          ELSE
     867     2104978 :             v = v + helium_spline(endpspline, arp)
     868             :          END IF
     869     6754137 :          IF (ar > 0.5_dp*helium%cell_size) THEN
     870     4649772 :             IF (nocut) v = v + helium_spline(endpspline, 0.5_dp*helium%cell_size)
     871             :          ELSE
     872     2104365 :             v = v + helium_spline(endpspline, ar)
     873             :          END IF
     874     6754137 :          res = 0.5_dp*v
     875             :       ELSE
     876             :          ! end-point action (first term):
     877    66763260 :          v = 0.5_dp*(helium_spline(endpspline, ar) + helium_spline(endpspline, arp))
     878   267053040 :          s = helium_pbc_r2(helium, r - rp)
     879    66763260 :          q = 0.5_dp*(ar + arp)
     880    66763260 :          z = (ar - arp)**2
     881    66763260 :          nsp = ((helium%pdx + 2)*(helium%pdx + 1))/2 - 1
     882    66763260 :          a = endpspline%x1
     883    66763260 :          b = endpspline%invh
     884    66763260 :          IF (q < a) THEN
     885           0 :             b = b*(q - a)
     886           0 :             a = 1.0_dp - b
     887             :             work(1:nsp) = a*offdiagsplines(1:nsp, 1, 1) + b*(offdiagsplines(1:nsp, 1, 2) - &
     888           0 :                                                              offdiagsplines(1:nsp, 2, 2)*endpspline%h26)
     889    66763260 :          ELSE IF (q > endpspline%xn) THEN
     890      764446 :             b = b*(q - endpspline%xn) + 1.0_dp
     891      764446 :             a = 1.0_dp - b
     892      764446 :             tline = endpspline%n
     893             :             work(1:nsp) = b*offdiagsplines(1:nsp, 1, tline) + a*(offdiagsplines(1:nsp, 1, tline - 1) - &
     894     4586676 :                                                                  offdiagsplines(1:nsp, 2, tline - 1)*endpspline%h26)
     895             :          ELSE
     896    65998814 :             a = (a - q)*b
     897    65998814 :             tline = INT(1.0_dp - a)
     898    65998814 :             a = a + REAL(tline, kind=dp)
     899    65998814 :             b = 1.0_dp - a
     900    65998814 :             h26 = -a*b*endpspline%h26
     901             :             work(1:nsp) = a*offdiagsplines(1:nsp, 1, tline) + b*offdiagsplines(1:nsp, 1, tline + 1) + &
     902             :                           h26*((a + 1.0_dp)*offdiagsplines(1:nsp, 2, tline) + (b + 1.0_dp)* &
     903   395992884 :                                offdiagsplines(1:nsp, 2, tline + 1))
     904             :          END IF
     905    66763260 :          work(nsp + 1) = v
     906    66763260 :          tline = 1
     907    66763260 :          v = work(1)
     908   200289780 :          DO i = 1, helium%pdx
     909   133526520 :             v = v*z
     910   133526520 :             tline = tline + 1
     911   133526520 :             b = work(tline)
     912   333816300 :             DO j = 1, i
     913   200289780 :                tline = tline + 1
     914   333816300 :                b = b*s + work(tline)
     915             :             END DO
     916   200289780 :             v = v + b
     917             :          END DO
     918             :          res = v
     919             :       END IF
     920    73517397 :    END FUNCTION helium_eval_expansion
     921             : 
     922             : ! ***************************************************************************
     923             : !> \brief  Calculate the pair-product action or energy by evaluating the
     924             : !>         power series expansion according to Eq. 4.46 in Ceperley 1995.
     925             : !> \param helium ...
     926             : !> \param rchain ...
     927             : !> \param nchain ...
     928             : !> \param aij ...
     929             : !> \param vcoef ...
     930             : !> \param energy ...
     931             : !> \return ...
     932             : !> \author Harald Forbert
     933             : !> \note This version calculates the action/energy for a chain segment
     934             : ! **************************************************************************************************
     935     5057417 :    FUNCTION helium_eval_chain(helium, rchain, nchain, aij, vcoef, energy) RESULT(res)
     936             : 
     937             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     938             :       INTEGER, INTENT(IN)                                :: nchain
     939             :       REAL(KIND=dp), DIMENSION(3, nchain), INTENT(INOUT) :: rchain
     940             :       REAL(KIND=dp), DIMENSION(nchain), INTENT(INOUT)    :: aij
     941             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: vcoef
     942             :       LOGICAL, INTENT(IN), OPTIONAL                      :: energy
     943             :       REAL(KIND=dp)                                      :: res
     944             : 
     945             :       INTEGER                                            :: chainpos, i, j, ndim, nsp, tline
     946             :       LOGICAL                                            :: nocut
     947             :       REAL(KIND=dp)                                      :: a, ar, arp, b, ch, h26, q, s, totalv, v, &
     948             :                                                             z
     949             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     950     5057417 :          POINTER                                         :: offdiagsplines
     951             :       TYPE(spline_data_type), POINTER                    :: endpspline
     952             : 
     953     5057417 :       endpspline => helium%u0
     954     5057417 :       offdiagsplines => helium%uoffdiag
     955     5057417 :       nocut = .TRUE.
     956     5057417 :       IF (PRESENT(energy)) THEN
     957     3289662 :          IF (energy) THEN
     958     3289662 :             endpspline => helium%e0
     959     3289662 :             offdiagsplines => helium%eoffdiag
     960     3289662 :             nocut = .FALSE.
     961             :          END IF
     962             :       END IF
     963             : 
     964     5057417 :       ndim = helium%pdx
     965     5057417 :       nsp = ((ndim + 2)*(ndim + 1))/2 - 1
     966    35401919 :       vcoef(1:nsp + 1) = 0.0_dp
     967    84366303 :       totalv = 0.0_dp
     968    84366303 :       DO i = 1, nchain
     969    84366303 :          aij(i) = SQRT(helium_pbc_r2(helium, rchain(:, i)))
     970             :       END DO
     971    79308886 :       DO i = 2, nchain
     972   297005876 :          rchain(:, i - 1) = rchain(:, i - 1) - rchain(:, i)
     973    79308886 :          rchain(1, i - 1) = helium_pbc_r2(helium, rchain(:, i - 1))
     974             :       END DO
     975             :       !
     976     5057417 :       IF (helium%periodic) THEN
     977     3319852 :          ch = 0.5_dp*helium%cell_size
     978     3319852 :          IF (nocut) THEN
     979      159247 :             DO i = 2, nchain - 1
     980      159247 :                totalv = totalv + helium_spline(endpspline, MIN(aij(i), ch))
     981             :             END DO
     982       41292 :             totalv = totalv + 0.5_dp*helium_spline(endpspline, MIN(aij(1), ch))
     983       41292 :             totalv = totalv + 0.5_dp*helium_spline(endpspline, MIN(aij(nchain), ch))
     984             :          ELSE
     985    67679200 :             DO i = 2, nchain - 1
     986    64400640 :                ar = aij(i)
     987    67679200 :                IF (ar < ch) totalv = totalv + helium_spline(endpspline, ar)
     988             :             END DO
     989     3278560 :             ar = aij(1)
     990     3278560 :             IF (ar < ch) totalv = totalv + 0.5_dp*helium_spline(endpspline, ar)
     991     3278560 :             ar = aij(nchain)
     992     3278560 :             IF (ar < ch) totalv = totalv + 0.5_dp*helium_spline(endpspline, ar)
     993             :          END IF
     994             :       ELSE
     995     6413022 :          DO i = 2, nchain - 1
     996     6413022 :             totalv = totalv + helium_spline(endpspline, aij(i))
     997             :          END DO
     998     1737565 :          totalv = totalv + 0.5_dp*helium_spline(endpspline, aij(1))
     999     1737565 :          totalv = totalv + 0.5_dp*helium_spline(endpspline, aij(nchain))
    1000             :       END IF
    1001             : 
    1002    79308886 :       DO chainpos = 1, nchain - 1
    1003    74251469 :          ar = aij(chainpos)
    1004    74251469 :          arp = aij(chainpos + 1)
    1005    74251469 :          IF (helium%periodic .AND. ((ar > 0.5_dp*helium%cell_size) &
    1006     5057417 :                                     .OR. (arp > 0.5_dp*helium%cell_size))) THEN
    1007             :             CYCLE
    1008             :          ELSE
    1009    67951527 :             q = 0.5_dp*(ar + arp)
    1010    67951527 :             s = rchain(1, chainpos)
    1011    67951527 :             z = (ar - arp)**2
    1012    67951527 :             a = endpspline%x1
    1013    67951527 :             b = endpspline%invh
    1014    67951527 :             IF (q < a) THEN
    1015         678 :                b = b*(q - a)
    1016         678 :                a = 1.0_dp - b
    1017             :                vcoef(1:nsp) = a*offdiagsplines(1:nsp, 1, 1) + b*(offdiagsplines(1:nsp, 1, 2) - &
    1018        4068 :                                                                  offdiagsplines(1:nsp, 2, 2)*endpspline%h26)
    1019    67950849 :             ELSE IF (q > endpspline%xn) THEN
    1020     5555970 :                b = b*(q - endpspline%xn) + 1.0_dp
    1021     5555970 :                a = 1.0_dp - b
    1022     5555970 :                tline = endpspline%n
    1023             :                vcoef(1:nsp) = b*offdiagsplines(1:nsp, 1, tline) + a*(offdiagsplines(1:nsp, 1, tline - 1) - &
    1024    33335820 :                                                                      offdiagsplines(1:nsp, 2, tline - 1)*endpspline%h26)
    1025             :             ELSE
    1026    62394879 :                a = (a - q)*b
    1027    62394879 :                tline = INT(1.0_dp - a)
    1028    62394879 :                a = a + REAL(tline, kind=dp)
    1029    62394879 :                b = 1.0_dp - a
    1030    62394879 :                h26 = -a*b*endpspline%h26
    1031             :                vcoef(1:nsp) = a*offdiagsplines(1:nsp, 1, tline) + b*offdiagsplines(1:nsp, 1, tline + 1) + &
    1032             :                               h26*((a + 1.0_dp)*offdiagsplines(1:nsp, 2, tline) + (b + 1.0_dp)* &
    1033   374369274 :                                    offdiagsplines(1:nsp, 2, tline + 1))
    1034             :             END IF
    1035    67951527 :             tline = 1
    1036    67951527 :             v = vcoef(1)
    1037   203854581 :             DO i = 1, ndim
    1038   135903054 :                tline = tline + 1
    1039   135903054 :                b = vcoef(tline)
    1040   339757635 :                DO j = 1, i
    1041   203854581 :                   tline = tline + 1
    1042   339757635 :                   b = b*s + vcoef(tline)
    1043             :                END DO
    1044   203854581 :                v = v*z + b
    1045             :             END DO
    1046    67951527 :             totalv = totalv + v
    1047             :          END IF
    1048             :       END DO
    1049     5057417 :       res = totalv
    1050     5057417 :    END FUNCTION helium_eval_chain
    1051             : 
    1052             : ! **************************************************************************************************
    1053             : !> \brief ...
    1054             : !> \param helium ...
    1055             : !> \author Harald Forbert
    1056             : ! **************************************************************************************************
    1057        6600 :    SUBROUTINE helium_update_transition_matrix(helium)
    1058             : 
    1059             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1060             : 
    1061             :       INTEGER                                            :: b, c, i, j, k, m, n, nb
    1062        6600 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: lens, order
    1063        6600 :       INTEGER, DIMENSION(:), POINTER                     :: perm
    1064        6600 :       INTEGER, DIMENSION(:, :), POINTER                  :: nmatrix
    1065             :       REAL(KIND=dp)                                      :: f, q, t, v
    1066        6600 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: p
    1067             :       REAL(KIND=dp), DIMENSION(3)                        :: r
    1068        6600 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ipmatrix, pmatrix, tmatrix
    1069        6600 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos
    1070             : 
    1071        6600 :       nb = helium%atoms
    1072             :       !TODO: check allocation status
    1073       19800 :       ALLOCATE (p(2*nb))
    1074       19800 :       ALLOCATE (order(nb))
    1075       19800 :       ALLOCATE (lens(2*nb))
    1076        6600 :       b = helium%beads - helium%bisection + 1
    1077        6600 :       f = -0.5_dp/(helium%hb2m*helium%tau*helium%bisection)
    1078        6600 :       tmatrix => helium%tmatrix
    1079        6600 :       pmatrix => helium%pmatrix
    1080        6600 :       ipmatrix => helium%ipmatrix
    1081        6600 :       nmatrix => helium%nmatrix
    1082        6600 :       perm => helium%permutation
    1083        6600 :       pos => helium%pos
    1084      217800 :       DO i = 1, nb
    1085     6969600 :          DO j = 1, nb
    1086     6758400 :             v = 0.0_dp
    1087    27033600 :             r(:) = pos(:, i, b) - pos(:, j, 1)
    1088     6758400 :             CALL helium_pbc(helium, r)
    1089     6758400 :             v = v + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
    1090     6969600 :             pmatrix(i, j) = f*v
    1091             :          END DO
    1092      211200 :          t = pmatrix(i, perm(i)) ! just some reference
    1093      211200 :          v = 0.0_dp
    1094     6969600 :          DO j = 1, nb
    1095     6758400 :             tmatrix(i, j) = EXP(pmatrix(i, j) - t)
    1096     6969600 :             v = v + tmatrix(i, j)
    1097             :          END DO
    1098             :          ! normalize
    1099             :          q = t + LOG(v)
    1100      211200 :          t = 1.0_dp/v
    1101     6969600 :          DO j = 1, nb
    1102     6758400 :             tmatrix(i, j) = tmatrix(i, j)*t
    1103     6969600 :             ipmatrix(i, j) = 1.0_dp/tmatrix(i, j)
    1104             :          END DO
    1105             : 
    1106             :          ! at this point we have:
    1107             :          ! tmatrix(i,j) = exp(-f*(r_i^b - r_j^1)**2) normalized such
    1108             :          !    that sum_j tmatrix(i,j) = 1.
    1109             :          !    ( tmatrix(k1,k2) = t_{k1,k2} / h_{k1} of ceperly. )
    1110             :          !    so tmatrix(i,j) is the probability to try to change a permutation
    1111             :          !    with particle j (assuming particle i is already selected as well)
    1112             :          ! ipmatrix(i,j) = 1.0/tmatrix(i,j)
    1113             :          ! pmatrix(i,j) = log(tmatrix(i,j))  + some_offset(i)
    1114             : 
    1115             :          ! generate optimal search tree so we can select which particle j
    1116             :          ! belongs to a given random_number as fast as possible.
    1117             :          ! (the traditional approach would be to generate a table
    1118             :          !  of cumulative probabilities and to search that table)
    1119             :          ! so for example if we have:
    1120             :          ! tmatrix(i,:) = ( 0.1 , 0.4 , 0.2 , 0.3 )
    1121             :          ! traditionally we would build the running sum table:
    1122             :          !  ( 0.1 , 0.5 , 0.7 , 1.0 ) and for a random number r
    1123             :          ! would search this table for the lowest index larger than r
    1124             :          ! (which would then be the particle index chosen by this random number)
    1125             :          ! we build an optimal binary search tree instead, so here
    1126             :          ! we would have:
    1127             :          ! if ( r > 0.6 ) then take index 2,
    1128             :          ! else if ( r > 0.3 ) then take index 4,
    1129             :          ! else if ( r > 0.1 ) then take index 3 else index 1.
    1130             :          ! the search tree is generated in tmatrix and nmatrix.
    1131             :          ! tmatrix contains the decision values (0.6,0.3,0.1 in this case)
    1132             :          ! and nmatrix contains the two branches (what to do if lower or higher)
    1133             :          ! negative numbers in nmatrix mean take minus that index
    1134             :          ! positive number means go down the tree to that next node, since we
    1135             :          ! put the root of the tree at the end the arrays in the example would
    1136             :          ! look like this:
    1137             :          ! tmatrix(i,:) = ( 0.1 , 0.3 , 0.6 , arbitrary )
    1138             :          ! namtrix(i,:) = ( -1 , -3 , 1 , -4 , 2 , -2 , arb. , arb. )
    1139             :          !
    1140             :          ! the way to generate this tree may not be the best, but the
    1141             :          ! tree generation itself shouldn't be needed quite that often:
    1142             :          !
    1143             :          ! first sort values (with some variant of heap sort)
    1144             : 
    1145     6969600 :          DO j = 1, nb
    1146     6758400 :             order(j) = j
    1147     6969600 :             p(j) = tmatrix(i, j)
    1148             :          END DO
    1149      211200 :          IF (nb > 1) THEN ! if nb = 1 it is already sorted.
    1150      211200 :             k = nb/2 + 1
    1151      211200 :             c = nb
    1152     9715200 :             DO
    1153     9926400 :                IF (k > 1) THEN
    1154             :                   ! building up the heap:
    1155     3379200 :                   k = k - 1
    1156     3379200 :                   n = order(k)
    1157     3379200 :                   v = p(k)
    1158             :                ELSE
    1159             :                   ! removing the top of the heap
    1160     6547200 :                   n = order(c)
    1161     6547200 :                   v = p(c)
    1162     6547200 :                   order(c) = order(1)
    1163     6547200 :                   p(c) = p(1)
    1164     6547200 :                   c = c - 1
    1165     6547200 :                   IF (c == 1) THEN
    1166      211200 :                      order(1) = n
    1167      211200 :                      p(1) = v
    1168             :                      EXIT
    1169             :                   END IF
    1170             :                END IF
    1171     9715200 :                m = k
    1172     9715200 :                j = 2*k
    1173             :                ! restoring heap order between k and c
    1174    21788037 :                DO
    1175    31503237 :                   IF (j > c) EXIT
    1176    24171708 :                   IF (j < c) THEN
    1177    22911836 :                      IF (p(j) < p(j + 1)) j = j + 1
    1178             :                   END IF
    1179    24171708 :                   IF (v >= p(j)) EXIT
    1180    21788037 :                   order(m) = order(j)
    1181    21788037 :                   p(m) = p(j)
    1182    21788037 :                   m = j
    1183    24171708 :                   j = 2*j
    1184             :                END DO
    1185     9715200 :                order(m) = n
    1186     9715200 :                p(m) = v
    1187             :             END DO
    1188             :          END IF
    1189             : 
    1190             :          ! now:
    1191             :          !    p(1:nb)     : tmatrix(i,1:nb) sorted in ascending order
    1192             :          !    order(1:nb) : corresponding index: p(j) == tmatrix(i,order(j))
    1193             :          !                                                       for all j
    1194             : 
    1195             :          ! merge sort with elements as we generate new interior search nodes
    1196             :          ! by combining older elements/nodes
    1197             : 
    1198             :          ! first fill unused part of array with guard values:
    1199     6969600 :          DO j = nb + 1, 2*nb
    1200     6969600 :             p(j) = 2.0_dp
    1201             :          END DO
    1202             : 
    1203             :          ! j   - head of leaf queue
    1204             :          ! c+1 - head of node queue in p (c in lens)
    1205             :          ! m+1 - tail of node queue in p (m in lens)
    1206             :          c = nb + 1
    1207             :          j = 1
    1208     6758400 :          DO m = nb + 1, 2*nb - 1
    1209             :             ! get next smallest element
    1210     6547200 :             IF (p(j) < p(c + 1)) THEN
    1211     2319018 :                v = p(j)
    1212     2319018 :                lens(j) = m
    1213     2319018 :                j = j + 1
    1214             :             ELSE
    1215     4228182 :                v = p(c + 1)
    1216     4228182 :                lens(c) = m
    1217     4228182 :                c = c + 1
    1218             :             END IF
    1219             :             ! get the second next smallest element
    1220     6758400 :             IF (p(j) < p(c + 1)) THEN
    1221     4439382 :                p(m + 1) = v + p(j)
    1222     4439382 :                lens(j) = m
    1223     4439382 :                j = j + 1
    1224             :             ELSE
    1225     2107818 :                p(m + 1) = v + p(c + 1)
    1226     2107818 :                lens(c) = m
    1227     2107818 :                c = c + 1
    1228             :             END IF
    1229             :          END DO
    1230             : 
    1231             :          ! lens(:) now has the tree with lens(j) pointing to its parent
    1232             :          ! the root of the tree is at 2*nb-1
    1233             :          ! calculate the depth of each node in the tree now: (root = 0)
    1234             : 
    1235      211200 :          lens(2*nb - 1) = 0
    1236    13305600 :          DO m = 2*nb - 2, 1, -1
    1237    13305600 :             lens(m) = lens(lens(m)) + 1
    1238             :          END DO
    1239             : 
    1240             :          ! lens(:) now has the depths of the nodes/leafs
    1241             : 
    1242             : #if 0
    1243             :          ! calculate average search depth (for information only)
    1244             :          v = 0.0_dp
    1245             :          DO j = 1, nb
    1246             :             v = v + p(j)*lens(j)
    1247             :          END DO
    1248             :          PRINT *, "Expected number of comparisons with i=", i, v
    1249             : #endif
    1250             : 
    1251             :          ! reset the nodes, for the canonical tree we just need the leaf info
    1252     6969600 :          DO j = 1, nb
    1253     6758400 :             lens(j + nb) = 0
    1254     6969600 :             p(j + nb) = 0.0_dp
    1255             :          END DO
    1256             : 
    1257             :          ! build the canonical tree (number of decisions on average are
    1258             :          ! the same to the tree we build above, but it has better caching behavior
    1259             : 
    1260             :          ! c head of leafs
    1261             :          ! m head of interior nodes
    1262             :          c = 1
    1263             :          m = nb + 1
    1264    13305600 :          DO k = 1, 2*nb - 2
    1265    13094400 :             j = nb + 1 + (k - 1)/2
    1266    13094400 :             IF (lens(c) > lens(m + 1)) THEN
    1267     6758400 :                nmatrix(i, k) = -order(c)
    1268     6758400 :                lens(j + 1) = lens(c) - 1
    1269     6758400 :                v = p(c)
    1270     6758400 :                c = c + 1
    1271             :             ELSE
    1272     6336000 :                nmatrix(i, k) = m - nb
    1273     6336000 :                lens(j + 1) = lens(m + 1) - 1
    1274     6336000 :                v = p(m)
    1275     6336000 :                m = m + 1
    1276             :             END IF
    1277    13094400 :             p(j) = p(j) + v
    1278    13305600 :             IF (MOD(k, 2) == 1) tmatrix(i, j - nb) = v
    1279             :          END DO
    1280             : 
    1281             :          ! now:
    1282             :          !    nmatrix(i,2*j+1) left child of node j
    1283             :          !    nmatrix(i,2*j+2) right child of node j
    1284             :          !       children:
    1285             :          !           negative : leaf with particle index == abs(value)
    1286             :          !           positive : child node index
    1287             :          !    p(j) weight of leaf j
    1288             :          !    p(nb+j) weight of node j
    1289             :          !    tmatrix(i,j) weight of left child of node j
    1290             : 
    1291             :          ! fix offsets for decision tree:
    1292             : 
    1293      211200 :          p(nb - 1) = 0.0_dp
    1294     6765000 :          DO m = nb - 1, 1, -1
    1295             :             ! if right child is a node, set its offset and
    1296             :             ! change its decision value
    1297     6547200 :             IF (nmatrix(i, 2*m) > 0) THEN
    1298      564170 :                p(nmatrix(i, 2*m)) = tmatrix(i, m)
    1299      564170 :                tmatrix(i, nmatrix(i, 2*m)) = tmatrix(i, nmatrix(i, 2*m)) + tmatrix(i, m)
    1300             :             END IF
    1301             :             ! if left child is a node, set its offset and
    1302             :             ! change its decision value
    1303     6758400 :             IF (nmatrix(i, 2*m - 1) > 0) THEN
    1304     5771830 :                p(nmatrix(i, 2*m - 1)) = p(m)
    1305     5771830 :                tmatrix(i, nmatrix(i, 2*m - 1)) = tmatrix(i, nmatrix(i, 2*m - 1)) + p(m)
    1306             :             END IF
    1307             :          END DO
    1308             : 
    1309             :          ! canonical optimal search tree done
    1310             : 
    1311             : #if 0
    1312             :          !some test code, to check if it gives the right distribution
    1313             :          DO k = 1, nb
    1314             :             p(k) = 1.0/ipmatrix(i, k)
    1315             :          END DO
    1316             :          lens(:) = 0
    1317             :          ! number of random numbers to generate:
    1318             :          c = 1000000000
    1319             :          DO j = 1, c
    1320             :             v = helium%rng_stream_uniform%next()
    1321             :             ! walk down the search tree:
    1322             :             k = nb - 1
    1323             :             DO
    1324             :                IF (tmatrix(i, k) > v) THEN
    1325             :                   k = nmatrix(i, 2*k - 1)
    1326             :                ELSE
    1327             :                   k = nmatrix(i, 2*k)
    1328             :                END IF
    1329             :                IF (k < 0) EXIT
    1330             :             END DO
    1331             :             k = -k
    1332             :             ! increment the counter for this particle index
    1333             :             lens(k) = lens(k) + 1
    1334             :          END DO
    1335             :          ! search for maximum deviation from expectation value
    1336             :          ! (relative to the expected variance)
    1337             :          v = 0.0_dp
    1338             :          k = -1
    1339             :          DO j = 1, nb
    1340             :             q = ABS((lens(j) - c*p(j))/SQRT(c*p(j)))
    1341             :             !PRINT *,j,lens(j),c*p(j)
    1342             :             IF (q > v) THEN
    1343             :                v = q
    1344             :                k = j
    1345             :             END IF
    1346             :             !PRINT *,lens(j),c*p(j),(lens(j)-c*p(j))/sqrt(c*p(j))
    1347             :          END DO
    1348             :          PRINT *, "MAXDEV:", k, lens(k), c*p(k), v
    1349             :          !PRINT *,"TMAT:",tmatrix(i,:)
    1350             :          !PRINT *,"NMAT:",nmatrix(i,:)
    1351             :          !STOP
    1352             : #endif
    1353             : #if 0
    1354             :          !additional test code:
    1355             :          p(:) = -1.0_dp
    1356             :          p(nb - 1) = 0.0_dp
    1357             :          p(2*nb - 1) = 1.0_dp
    1358             :          DO j = nb - 1, 1, -1
    1359             :             ! right child
    1360             :             IF (nmatrix(i, 2*j) > 0) THEN
    1361             :                c = nmatrix(i, 2*j)
    1362             :                p(c) = tmatrix(i, j)
    1363             :                p(c + nb) = p(j + nb)
    1364             :             ELSE
    1365             :                c = -nmatrix(i, 2*j)
    1366             :                !PRINT *,c,1.0/ipmatrix(i,c),p(j+nb)-tmatrix(i,j)
    1367             :                IF (ABS(1.0/ipmatrix(i, c) - (p(j + nb) - tmatrix(i, j))) > &
    1368             :                    10.0_dp*EPSILON(1.0_dp)) THEN
    1369             :                   PRINT *, "Probability mismatch for particle i->j", i, c
    1370             :                   PRINT *, "Got", p(j + nb) - tmatrix(i, j), "should be", 1.0/ipmatrix(i, c)
    1371             :                   STOP
    1372             :                END IF
    1373             :             END IF
    1374             :             ! left child
    1375             :             IF (nmatrix(i, 2*j - 1) > 0) THEN
    1376             :                c = nmatrix(i, 2*j - 1)
    1377             :                p(c + nb) = tmatrix(i, j)
    1378             :                p(c) = p(j)
    1379             :             ELSE
    1380             :                c = -nmatrix(i, 2*j - 1)
    1381             :                !PRINT *,c,1.0/ipmatrix(i,c),tmatrix(i,j)-p(j)
    1382             :                IF (ABS(1.0/ipmatrix(i, c) - (tmatrix(i, j) - p(j))) > &
    1383             :                    10.0_dp*EPSILON(1.0_dp)) THEN
    1384             :                   PRINT *, "Probability mismatch for particle i->j", i, c
    1385             :                   PRINT *, "Got", tmatrix(i, j) - p(j), "should be", 1.0/ipmatrix(i, c)
    1386             :                   STOP
    1387             :                END IF
    1388             :             END IF
    1389             :          END DO
    1390             :          PRINT *, "Probabilities ok"
    1391             : #endif
    1392             : 
    1393             :       END DO
    1394             : 
    1395             :       ! initialize trial permutation with some identity permutation
    1396             :       ! (should not be taken, but just in case it does we have something valid)
    1397             : 
    1398        6600 :       helium%pweight = 0.0_dp
    1399        6600 :       t = helium%rng_stream_uniform%next()
    1400        6600 :       helium%ptable(1) = 1 + INT(t*nb)
    1401        6600 :       helium%ptable(2) = -1
    1402             : 
    1403             :       ! recalculate inverse permutation table (just in case)
    1404      217800 :       DO i = 1, nb
    1405      217800 :          helium%iperm(perm(i)) = i
    1406             :       END DO
    1407             : 
    1408             :       ! clean up:
    1409        6600 :       DEALLOCATE (lens)
    1410        6600 :       DEALLOCATE (order)
    1411        6600 :       DEALLOCATE (p)
    1412             : 
    1413        6600 :    END SUBROUTINE helium_update_transition_matrix
    1414             : 
    1415             : ! **************************************************************************************************
    1416             : !> \brief ...
    1417             : !> \param spl ...
    1418             : !> \param xx ...
    1419             : !> \return ...
    1420             : !> \author Harald Forbert
    1421             : ! **************************************************************************************************
    1422   342551603 :    FUNCTION helium_spline(spl, xx) RESULT(res)
    1423             :       TYPE(spline_data_type), INTENT(IN)                 :: spl
    1424             :       REAL(KIND=dp), INTENT(IN)                          :: xx
    1425             :       REAL(KIND=dp)                                      :: res
    1426             : 
    1427             :       REAL(KIND=dp)                                      :: a, b
    1428             : 
    1429   342551603 :       IF (xx < spl%x1) THEN
    1430        1995 :          b = spl%invh*(xx - spl%x1)
    1431        1995 :          a = 1.0_dp - b
    1432        1995 :          res = a*spl%y(1) + b*(spl%y(2) - spl%y2(2)*spl%h26)
    1433   342549608 :       ELSE IF (xx > spl%xn) THEN
    1434     8988646 :          b = spl%invh*(xx - spl%xn) + 1.0_dp
    1435     8988646 :          a = 1.0_dp - b
    1436     8988646 :          res = b*spl%y(spl%n) + a*(spl%y(spl%n - 1) - spl%y2(spl%n - 1)*spl%h26)
    1437             :       ELSE
    1438   333560962 :          res = spline_value(spl, xx)
    1439             :       END IF
    1440   342551603 :    END FUNCTION helium_spline
    1441             : 
    1442             : ! ***************************************************************************
    1443             : !> \brief  Return the distance <rij> between bead <ib> of atom <ia>
    1444             : !>         and bead <jb> of atom <ja>.
    1445             : !> \param helium ...
    1446             : !> \param ia ...
    1447             : !> \param ib ...
    1448             : !> \param ja ...
    1449             : !> \param jb ...
    1450             : !> \return ...
    1451             : !> \date   2009-07-17
    1452             : !> \author Lukasz Walewski
    1453             : ! **************************************************************************************************
    1454           0 :    PURE FUNCTION helium_bead_rij(helium, ia, ib, ja, jb) RESULT(rij)
    1455             : 
    1456             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1457             :       INTEGER, INTENT(IN)                                :: ia, ib, ja, jb
    1458             :       REAL(KIND=dp)                                      :: rij
    1459             : 
    1460             :       REAL(KIND=dp)                                      :: dx, dy, dz
    1461             : 
    1462           0 :       dx = helium%pos(1, ia, ib) - helium%pos(1, ja, jb)
    1463           0 :       dy = helium%pos(2, ia, ib) - helium%pos(2, ja, jb)
    1464           0 :       dz = helium%pos(3, ia, ib) - helium%pos(3, ja, jb)
    1465           0 :       rij = SQRT(dx*dx + dy*dy + dz*dz)
    1466             : 
    1467           0 :    END FUNCTION helium_bead_rij
    1468             : 
    1469             : ! ***************************************************************************
    1470             : !> \brief  Given the atom number and permutation state return the cycle
    1471             : !>         number the atom belongs to.
    1472             : !> \param helium ...
    1473             : !> \param atom_number ...
    1474             : !> \param permutation ...
    1475             : !> \return ...
    1476             : !> \date   2009-07-21
    1477             : !> \author Lukasz Walewski
    1478             : !> \note   Cycles (or paths) are numbered from 1 to <num_cycles>, where
    1479             : !>         <num_cycles> is in the range of (1, <helium%atoms>).
    1480             : !>         if (num_cycles .EQ. 1) then all atoms belong to one cycle
    1481             : !>         if (num_cycles .EQ. helium%atoms) then there are no cycles of
    1482             : !>         length greater than 1 (i.e. no atoms are connected)
    1483             : ! **************************************************************************************************
    1484         250 :    FUNCTION helium_cycle_number(helium, atom_number, permutation) RESULT(cycle_number)
    1485             : 
    1486             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1487             :       INTEGER, INTENT(IN)                                :: atom_number
    1488             :       INTEGER, DIMENSION(:), POINTER                     :: permutation
    1489             :       INTEGER                                            :: cycle_number
    1490             : 
    1491             :       INTEGER                                            :: atom_idx, cycle_idx, cycle_num, ia, ib, &
    1492             :                                                             ic, num_cycles
    1493         250 :       INTEGER, DIMENSION(:), POINTER                     :: cycle_index
    1494             :       LOGICAL                                            :: found, new_cycle
    1495             : 
    1496         250 :       NULLIFY (cycle_index)
    1497         250 :       cycle_index => helium%itmp_atoms_1d
    1498        1500 :       cycle_index(:) = 0
    1499             : 
    1500         250 :       num_cycles = 0
    1501         250 :       found = .FALSE.
    1502         250 :       cycle_num = -1
    1503        1000 :       DO ia = 1, helium%atoms
    1504             :          ! this loop reaches its maximum iteration count when atom in question
    1505             :          ! is the last one (i.e. when atom_number .EQ. helium%atoms)
    1506             : 
    1507             :          ! exit if we have found the cycle number for the atom in question
    1508         950 :          IF (found) THEN
    1509             :             EXIT
    1510             :          END IF
    1511             : 
    1512             :          ! initialize current cycle index with the current atom
    1513         750 :          cycle_idx = ia
    1514             : 
    1515         750 :          atom_idx = ia
    1516        1000 :          DO ib = 1, helium%atoms*helium%beads
    1517             :             ! this loop reaches its maximum iteration count when all He atoms
    1518             :             ! form one cycle (i.e. all beads belong to one path)
    1519             : 
    1520             :             ! proceed along the path
    1521         750 :             atom_idx = permutation(atom_idx)
    1522             : 
    1523         750 :             IF (atom_idx .EQ. ia) THEN
    1524             :                ! end of cycle detected (looped back to the first atom)
    1525             : 
    1526             :                ! check if this is a new cycle
    1527             :                new_cycle = .TRUE.
    1528        1750 :                DO ic = 1, num_cycles
    1529        1750 :                   IF (cycle_index(ic) .EQ. cycle_idx) THEN
    1530           0 :                      new_cycle = .FALSE.
    1531             :                   END IF
    1532             :                END DO
    1533             : 
    1534         750 :                IF (new_cycle) THEN
    1535             :                   ! increase number of cycles and update the current cycle's index
    1536         750 :                   num_cycles = num_cycles + 1
    1537         750 :                   cycle_index(num_cycles) = cycle_idx
    1538             :                END IF
    1539             : 
    1540             :                ! if this was the atom in question
    1541         750 :                IF (ia .EQ. atom_number) THEN
    1542             :                   ! save the cycle index it belongs to
    1543         250 :                   cycle_num = cycle_idx
    1544             : 
    1545             :                   ! exit the loop over atoms, we've found what we've been looking for
    1546         250 :                   found = .TRUE.
    1547             :                END IF
    1548             : 
    1549             :                ! exit the loop over beads, there are no more (new) beads in this cycle
    1550             :                EXIT
    1551             :             END IF
    1552             : 
    1553             :             ! set the cycle index to the lowest atom index in this cycle
    1554           0 :             IF (atom_idx .LT. cycle_idx) THEN
    1555             :                cycle_idx = atom_idx
    1556             :             END IF
    1557             : 
    1558             :          END DO
    1559             : 
    1560             :       END DO
    1561             : 
    1562             : !TODO make it cp2k way
    1563         250 :       IF (.NOT. found) THEN
    1564           0 :          CPWARN("helium_cycle_number: we are going to return -1, problems ahead!")
    1565             :       END IF
    1566             : 
    1567             :       ! at this point we know the cycle index for atom <atom_number>
    1568             :       ! but it is expressed as the atom number of the first atom in that cycle
    1569             : 
    1570             :       ! renumber cycle indices, so that they form a range (1, <num_cycles>)
    1571             :       ! (don't do it actually - just return the requested <cycle_number>)
    1572         250 :       cycle_number = 0
    1573         750 :       DO ic = 1, num_cycles
    1574         750 :          IF (cycle_index(ic) .EQ. cycle_num) THEN
    1575             :             cycle_number = ic
    1576             :             EXIT
    1577             :          END IF
    1578             :       END DO
    1579             : 
    1580         250 :       NULLIFY (cycle_index)
    1581             : 
    1582         250 :    END FUNCTION helium_cycle_number
    1583             : 
    1584             : ! ***************************************************************************
    1585             : !> \brief  Given the atom number and permutation state return the length of
    1586             : !>         the path this atom belongs to.
    1587             : !> \param helium ...
    1588             : !> \param atom_number ...
    1589             : !> \param permutation ...
    1590             : !> \return ...
    1591             : !> \date   2009-10-07
    1592             : !> \author Lukasz Walewski
    1593             : ! **************************************************************************************************
    1594         250 :    PURE FUNCTION helium_path_length(helium, atom_number, permutation) RESULT(path_length)
    1595             : 
    1596             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1597             :       INTEGER, INTENT(IN)                                :: atom_number
    1598             :       INTEGER, DIMENSION(:), POINTER                     :: permutation
    1599             :       INTEGER                                            :: path_length
    1600             : 
    1601             :       INTEGER                                            :: atom_idx, ia
    1602             :       LOGICAL                                            :: path_end_reached
    1603             : 
    1604         250 :       atom_idx = atom_number
    1605         250 :       path_length = 0
    1606         250 :       path_end_reached = .FALSE.
    1607         250 :       DO ia = 1, helium%atoms
    1608         250 :          path_length = path_length + 1
    1609         250 :          atom_idx = permutation(atom_idx)
    1610         250 :          IF (atom_idx .EQ. atom_number) THEN
    1611             :             path_end_reached = .TRUE.
    1612             :             EXIT
    1613             :          END IF
    1614             :       END DO
    1615             : 
    1616         250 :       IF (.NOT. path_end_reached) THEN
    1617           0 :          path_length = -1
    1618             :       END IF
    1619             : 
    1620         250 :    END FUNCTION helium_path_length
    1621             : 
    1622             : ! ***************************************************************************
    1623             : !> \brief  Given an element of a permutation return the cycle it belongs to.
    1624             : !> \param element ...
    1625             : !> \param permutation ...
    1626             : !> \return ...
    1627             : !> \date   2013-12-10
    1628             : !> \author Lukasz Walewski
    1629             : !> \note   This function allocates memory and returns a pointer to it,
    1630             : !>         do not forget to deallocate once finished with the results.
    1631             : ! **************************************************************************************************
    1632         250 :    PURE FUNCTION helium_cycle_of(element, permutation) RESULT(CYCLE)
    1633             : 
    1634             :       INTEGER, INTENT(IN)                                :: element
    1635             :       INTEGER, DIMENSION(:), INTENT(IN), POINTER         :: permutation
    1636             :       INTEGER, DIMENSION(:), POINTER                     :: CYCLE
    1637             : 
    1638             :       INTEGER                                            :: ia, icur, len, nsize
    1639         250 :       INTEGER, DIMENSION(:), POINTER                     :: my_cycle
    1640             :       LOGICAL                                            :: cycle_end_reached
    1641             : 
    1642         250 :       nsize = SIZE(permutation)
    1643             : 
    1644             :       ! maximum possible cycle length is the number of atoms
    1645         250 :       NULLIFY (my_cycle)
    1646         750 :       ALLOCATE (my_cycle(nsize))
    1647             : 
    1648             :       ! traverse the permutation table
    1649         250 :       len = 0
    1650         250 :       icur = element
    1651         250 :       cycle_end_reached = .FALSE.
    1652         250 :       DO ia = 1, nsize
    1653         250 :          len = len + 1
    1654         250 :          my_cycle(len) = icur
    1655         250 :          icur = permutation(icur)
    1656         250 :          IF (icur .EQ. element) THEN
    1657             :             cycle_end_reached = .TRUE.
    1658             :             EXIT
    1659             :          END IF
    1660             :       END DO
    1661             : 
    1662         250 :       IF (.NOT. cycle_end_reached) THEN
    1663             :          ! return null
    1664           0 :          NULLIFY (CYCLE)
    1665             :       ELSE
    1666             :          ! assign the result
    1667         750 :          ALLOCATE (CYCLE(len))
    1668        1000 :          CYCLE(1:len) = my_cycle(1:len)
    1669             :       END IF
    1670             : 
    1671             :       ! clean up
    1672         250 :       DEALLOCATE (my_cycle)
    1673         250 :       NULLIFY (my_cycle)
    1674             : 
    1675         250 :    END FUNCTION helium_cycle_of
    1676             : 
    1677             : ! ***************************************************************************
    1678             : !> \brief  Return total winding number
    1679             : !> \param helium ...
    1680             : !> \return ...
    1681             : !> \date   2014-04-24
    1682             : !> \author Lukasz Walewski
    1683             : ! **************************************************************************************************
    1684         153 :    FUNCTION helium_total_winding_number(helium) RESULT(wnum)
    1685             : 
    1686             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1687             :       REAL(KIND=dp), DIMENSION(3)                        :: wnum
    1688             : 
    1689             :       INTEGER                                            :: ia, ib
    1690             :       REAL(KIND=dp), DIMENSION(3)                        :: rcur
    1691         153 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ri, rj
    1692             : 
    1693         612 :       wnum(:) = 0.0_dp
    1694        3369 :       DO ia = 1, helium%atoms
    1695             :          ! sum of contributions from the rest of bead pairs
    1696       64128 :          DO ib = 1, helium%beads - 1
    1697       60912 :             ri => helium%pos(:, ia, ib)
    1698       60912 :             rj => helium%pos(:, ia, ib + 1)
    1699      243648 :             rcur(:) = ri(:) - rj(:)
    1700       60912 :             CALL helium_pbc(helium, rcur)
    1701      246864 :             wnum(:) = wnum(:) + rcur(:)
    1702             :          END DO
    1703             :          ! contribution from the last and the first bead
    1704        3216 :          ri => helium%pos(:, ia, helium%beads)
    1705        3216 :          rj => helium%pos(:, helium%permutation(ia), 1)
    1706       12864 :          rcur(:) = ri(:) - rj(:)
    1707        3216 :          CALL helium_pbc(helium, rcur)
    1708       13017 :          wnum(:) = wnum(:) + rcur(:)
    1709             :       END DO
    1710             : 
    1711         153 :    END FUNCTION helium_total_winding_number
    1712             : 
    1713             : ! ***************************************************************************
    1714             : !> \brief  Return link winding number
    1715             : !> \param helium ...
    1716             : !> \param ia ...
    1717             : !> \param ib ...
    1718             : !> \return ...
    1719             : !> \date   2014-04-24
    1720             : !> \author Lukasz Walewski
    1721             : ! **************************************************************************************************
    1722           0 :    FUNCTION helium_link_winding_number(helium, ia, ib) RESULT(wnum)
    1723             : 
    1724             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1725             :       INTEGER, INTENT(IN)                                :: ia, ib
    1726             :       REAL(KIND=dp), DIMENSION(3)                        :: wnum
    1727             : 
    1728             :       INTEGER                                            :: ja1, ja2, jb1, jb2
    1729           0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: r1, r2
    1730             : 
    1731           0 :       IF (ib .EQ. helium%beads) THEN
    1732           0 :          ja1 = ia
    1733           0 :          ja2 = helium%permutation(ia)
    1734           0 :          jb1 = ib
    1735           0 :          jb2 = 1
    1736             :       ELSE
    1737           0 :          ja1 = ia
    1738           0 :          ja2 = ia
    1739           0 :          jb1 = ib
    1740           0 :          jb2 = ib + 1
    1741             :       END IF
    1742           0 :       r1 => helium%pos(:, ja1, jb1)
    1743           0 :       r2 => helium%pos(:, ja2, jb2)
    1744           0 :       wnum(:) = r1(:) - r2(:)
    1745           0 :       CALL helium_pbc(helium, wnum)
    1746             : 
    1747           0 :    END FUNCTION helium_link_winding_number
    1748             : 
    1749             : ! ***************************************************************************
    1750             : !> \brief  Return total winding number (sum over all links)
    1751             : !> \param helium ...
    1752             : !> \return ...
    1753             : !> \date   2014-04-24
    1754             : !> \author Lukasz Walewski
    1755             : !> \note   mostly for sanity checks
    1756             : ! **************************************************************************************************
    1757           0 :    FUNCTION helium_total_winding_number_linkwise(helium) RESULT(wnum)
    1758             : 
    1759             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1760             :       REAL(KIND=dp), DIMENSION(3)                        :: wnum
    1761             : 
    1762             :       INTEGER                                            :: ia, ib
    1763             : 
    1764           0 :       wnum(:) = 0.0_dp
    1765           0 :       DO ia = 1, helium%atoms
    1766           0 :          DO ib = 1, helium%beads
    1767           0 :             wnum(:) = wnum(:) + helium_link_winding_number(helium, ia, ib)
    1768             :          END DO
    1769             :       END DO
    1770             : 
    1771           0 :    END FUNCTION helium_total_winding_number_linkwise
    1772             : 
    1773             : ! ***************************************************************************
    1774             : !> \brief  Return cycle winding number
    1775             : !> \param helium ...
    1776             : !> \param CYCLE ...
    1777             : !> \param pos ...
    1778             : !> \return ...
    1779             : !> \date   2014-04-24
    1780             : !> \author Lukasz Walewski
    1781             : ! **************************************************************************************************
    1782         500 :    FUNCTION helium_cycle_winding_number(helium, CYCLE, pos) RESULT(wnum)
    1783             : 
    1784             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1785             :       INTEGER, DIMENSION(:), POINTER                     :: CYCLE
    1786             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos
    1787             :       REAL(KIND=dp), DIMENSION(3)                        :: wnum
    1788             : 
    1789             :       INTEGER                                            :: i1, i2, ia, ib, nsize
    1790             :       REAL(KIND=dp), DIMENSION(3)                        :: rcur
    1791         250 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ri, rj
    1792             : 
    1793         250 :       nsize = SIZE(CYCLE)
    1794             : 
    1795             :       ! traverse the path
    1796        1000 :       wnum(:) = 0.0_dp
    1797         500 :       DO ia = 1, nsize
    1798             :          ! contributions from all bead pairs of the current atom
    1799        4000 :          DO ib = 1, helium%beads - 1
    1800        3750 :             ri => pos(:, CYCLE(ia), ib)
    1801        3750 :             rj => pos(:, CYCLE(ia), ib + 1)
    1802       15000 :             rcur(:) = ri(:) - rj(:)
    1803        3750 :             CALL helium_pbc(helium, rcur)
    1804       15250 :             wnum(:) = wnum(:) + rcur(:)
    1805             :          END DO
    1806             :          ! contribution from the last bead of the current atom
    1807             :          ! and the first bead of the next atom
    1808         250 :          i1 = CYCLE(ia)
    1809         250 :          IF (ia .EQ. nsize) THEN
    1810         250 :             i2 = CYCLE(1)
    1811             :          ELSE
    1812           0 :             i2 = CYCLE(ia + 1)
    1813             :          END IF
    1814         250 :          ri => pos(:, i1, helium%beads)
    1815         250 :          rj => pos(:, i2, 1)
    1816        1000 :          rcur(:) = ri(:) - rj(:)
    1817         250 :          CALL helium_pbc(helium, rcur)
    1818        1250 :          wnum(:) = wnum(:) + rcur(:)
    1819             :       END DO
    1820             : 
    1821         250 :    END FUNCTION helium_cycle_winding_number
    1822             : 
    1823             : ! ***************************************************************************
    1824             : !> \brief  Return total winding number (sum over all cycles)
    1825             : !> \param helium ...
    1826             : !> \return ...
    1827             : !> \date   2014-04-24
    1828             : !> \author Lukasz Walewski
    1829             : !> \note   mostly for sanity checks
    1830             : ! **************************************************************************************************
    1831           0 :    FUNCTION helium_total_winding_number_cyclewise(helium) RESULT(wnum)
    1832             : 
    1833             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1834             :       REAL(KIND=dp), DIMENSION(3)                        :: wnum
    1835             : 
    1836             :       INTEGER                                            :: ic
    1837             :       REAL(KIND=dp), DIMENSION(3)                        :: wn
    1838           0 :       TYPE(int_arr_ptr), DIMENSION(:), POINTER           :: cycles
    1839             : 
    1840             : ! decompose the current permutation state into permutation cycles
    1841             : 
    1842           0 :       NULLIFY (cycles)
    1843           0 :       cycles => helium_calc_cycles(helium%permutation)
    1844             : 
    1845           0 :       wnum(:) = 0.0_dp
    1846           0 :       DO ic = 1, SIZE(cycles)
    1847           0 :          wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
    1848           0 :          wnum(:) = wnum(:) + wn(:)
    1849             :       END DO
    1850             : 
    1851           0 :       DEALLOCATE (cycles)
    1852             : 
    1853           0 :    END FUNCTION helium_total_winding_number_cyclewise
    1854             : 
    1855             : ! ***************************************************************************
    1856             : !> \brief  Return total projected area
    1857             : !> \param helium ...
    1858             : !> \return ...
    1859             : !> \date   2014-04-24
    1860             : !> \author Lukasz Walewski
    1861             : ! **************************************************************************************************
    1862         153 :    FUNCTION helium_total_projected_area(helium) RESULT(area)
    1863             : 
    1864             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1865             :       REAL(KIND=dp), DIMENSION(3)                        :: area
    1866             : 
    1867             :       INTEGER                                            :: ia, ib
    1868             :       REAL(KIND=dp), DIMENSION(3)                        :: r1, r12, r2, rcur
    1869             : 
    1870         612 :       area(:) = 0.0_dp
    1871        3369 :       DO ia = 1, helium%atoms
    1872             :          ! contributions from all links of the current atom
    1873       64128 :          DO ib = 1, helium%beads - 1
    1874      243648 :             r1(:) = helium%pos(:, ia, ib)
    1875      243648 :             r2(:) = helium%pos(:, ia, ib + 1)
    1876             :             ! comment out for non-PBC version -->
    1877      243648 :             r12(:) = r2(:) - r1(:)
    1878       60912 :             CALL helium_pbc(helium, r1)
    1879       60912 :             CALL helium_pbc(helium, r12)
    1880      243648 :             r2(:) = r1(:) + r12(:)
    1881             :             ! comment out for non-PBC version <--
    1882       60912 :             rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
    1883       60912 :             rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
    1884       60912 :             rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
    1885      246864 :             area(:) = area(:) + rcur(:)
    1886             :          END DO
    1887             :          ! contribution from the last bead of the current atom
    1888             :          ! and the first bead of the next atom
    1889       12864 :          r1(:) = helium%pos(:, ia, helium%beads)
    1890       12864 :          r2(:) = helium%pos(:, helium%permutation(ia), 1)
    1891             :          ! comment out for non-PBC version -->
    1892       12864 :          r12(:) = r2(:) - r1(:)
    1893        3216 :          CALL helium_pbc(helium, r1)
    1894        3216 :          CALL helium_pbc(helium, r12)
    1895       12864 :          r2(:) = r1(:) + r12(:)
    1896             :          ! comment out for non-PBC version <--
    1897        3216 :          rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
    1898        3216 :          rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
    1899        3216 :          rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
    1900       13017 :          area(:) = area(:) + rcur(:)
    1901             :       END DO
    1902         612 :       area(:) = 0.5_dp*area(:)
    1903             : 
    1904         153 :    END FUNCTION helium_total_projected_area
    1905             : 
    1906             : ! ***************************************************************************
    1907             : !> \brief  Return link projected area
    1908             : !> \param helium ...
    1909             : !> \param ia ...
    1910             : !> \param ib ...
    1911             : !> \return ...
    1912             : !> \date   2014-04-24
    1913             : !> \author Lukasz Walewski
    1914             : ! **************************************************************************************************
    1915           0 :    FUNCTION helium_link_projected_area(helium, ia, ib) RESULT(area)
    1916             : 
    1917             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1918             :       INTEGER                                            :: ia, ib
    1919             :       REAL(KIND=dp), DIMENSION(3)                        :: area
    1920             : 
    1921             :       INTEGER                                            :: ja1, ja2, jb1, jb2
    1922             :       REAL(KIND=dp), DIMENSION(3)                        :: r1, r12, r2
    1923             : 
    1924           0 :       IF (ib .EQ. helium%beads) THEN
    1925           0 :          ja1 = ia
    1926           0 :          ja2 = helium%permutation(ia)
    1927           0 :          jb1 = ib
    1928           0 :          jb2 = 1
    1929             :       ELSE
    1930           0 :          ja1 = ia
    1931           0 :          ja2 = ia
    1932           0 :          jb1 = ib
    1933           0 :          jb2 = ib + 1
    1934             :       END IF
    1935           0 :       r1(:) = helium%pos(:, ja1, jb1)
    1936           0 :       r2(:) = helium%pos(:, ja2, jb2)
    1937             :       ! comment out for non-PBC version -->
    1938           0 :       r12(:) = r2(:) - r1(:)
    1939           0 :       CALL helium_pbc(helium, r1)
    1940           0 :       CALL helium_pbc(helium, r12)
    1941           0 :       r2(:) = r1(:) + r12(:)
    1942             :       ! comment out for non-PBC version <--
    1943           0 :       area(1) = r1(2)*r2(3) - r1(3)*r2(2)
    1944           0 :       area(2) = r1(3)*r2(1) - r1(1)*r2(3)
    1945           0 :       area(3) = r1(1)*r2(2) - r1(2)*r2(1)
    1946           0 :       area(:) = 0.5_dp*area(:)
    1947             : 
    1948           0 :    END FUNCTION helium_link_projected_area
    1949             : 
    1950             : ! ***************************************************************************
    1951             : !> \brief  Return total projected area (sum over all links)
    1952             : !> \param helium ...
    1953             : !> \return ...
    1954             : !> \date   2014-04-24
    1955             : !> \author Lukasz Walewski
    1956             : !> \note   mostly for sanity checks
    1957             : ! **************************************************************************************************
    1958           0 :    FUNCTION helium_total_projected_area_linkwise(helium) RESULT(area)
    1959             : 
    1960             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1961             :       REAL(KIND=dp), DIMENSION(3)                        :: area
    1962             : 
    1963             :       INTEGER                                            :: ia, ib
    1964             : 
    1965           0 :       area(:) = 0.0_dp
    1966           0 :       DO ia = 1, helium%atoms
    1967           0 :          DO ib = 1, helium%beads
    1968           0 :             area(:) = area(:) + helium_link_projected_area(helium, ia, ib)
    1969             :          END DO
    1970             :       END DO
    1971             : 
    1972           0 :    END FUNCTION helium_total_projected_area_linkwise
    1973             : 
    1974             : ! ***************************************************************************
    1975             : !> \brief  Return cycle projected area
    1976             : !> \param helium ...
    1977             : !> \param CYCLE ...
    1978             : !> \return ...
    1979             : !> \date   2014-04-24
    1980             : !> \author Lukasz Walewski
    1981             : ! **************************************************************************************************
    1982           0 :    FUNCTION helium_cycle_projected_area(helium, CYCLE) RESULT(area)
    1983             : 
    1984             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    1985             :       INTEGER, DIMENSION(:), POINTER                     :: CYCLE
    1986             :       REAL(KIND=dp), DIMENSION(3)                        :: area
    1987             : 
    1988             :       INTEGER                                            :: i1, i2, ia, ib, nsize
    1989             :       REAL(KIND=dp), DIMENSION(3)                        :: rcur, rsum
    1990           0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ri, rj
    1991             : 
    1992           0 :       nsize = SIZE(CYCLE)
    1993             : 
    1994             :       ! traverse the path
    1995           0 :       rsum(:) = 0.0_dp
    1996           0 :       DO ia = 1, nsize
    1997             :          ! contributions from all bead pairs of the current atom
    1998           0 :          DO ib = 1, helium%beads - 1
    1999           0 :             ri => helium%pos(:, CYCLE(ia), ib)
    2000           0 :             rj => helium%pos(:, CYCLE(ia), ib + 1)
    2001           0 :             rcur(1) = ri(2)*rj(3) - ri(3)*rj(2)
    2002           0 :             rcur(2) = ri(3)*rj(1) - ri(1)*rj(3)
    2003           0 :             rcur(3) = ri(1)*rj(2) - ri(2)*rj(1)
    2004           0 :             rsum(:) = rsum(:) + rcur(:)
    2005             :          END DO
    2006             :          ! contribution from the last bead of the current atom
    2007             :          ! and the first bead of the next atom
    2008           0 :          i1 = CYCLE(ia)
    2009           0 :          IF (ia .EQ. nsize) THEN
    2010           0 :             i2 = CYCLE(1)
    2011             :          ELSE
    2012           0 :             i2 = CYCLE(ia + 1)
    2013             :          END IF
    2014           0 :          ri => helium%pos(:, i1, helium%beads)
    2015           0 :          rj => helium%pos(:, i2, 1)
    2016           0 :          rcur(1) = ri(2)*rj(3) - ri(3)*rj(2)
    2017           0 :          rcur(2) = ri(3)*rj(1) - ri(1)*rj(3)
    2018           0 :          rcur(3) = ri(1)*rj(2) - ri(2)*rj(1)
    2019           0 :          rsum(:) = rsum(:) + rcur(:)
    2020             :       END DO
    2021           0 :       area(:) = 0.5_dp*rsum(:)
    2022             : 
    2023           0 :    END FUNCTION helium_cycle_projected_area
    2024             : 
    2025             : ! ***************************************************************************
    2026             : !> \brief  Return cycle projected area (sum over all links)
    2027             : !> \param helium ...
    2028             : !> \param CYCLE ...
    2029             : !> \return ...
    2030             : !> \date   2014-04-24
    2031             : !> \author Lukasz Walewski
    2032             : !> \note   mostly for sanity checks
    2033             : ! **************************************************************************************************
    2034           0 :    FUNCTION helium_cycle_projected_area_pbc(helium, CYCLE) RESULT(area)
    2035             : 
    2036             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2037             :       INTEGER, DIMENSION(:), POINTER                     :: CYCLE
    2038             :       REAL(KIND=dp), DIMENSION(3)                        :: area
    2039             : 
    2040             :       INTEGER                                            :: i1, i2, ia, ib, nsize
    2041             :       REAL(KIND=dp), DIMENSION(3)                        :: r1, r12, r2, rcur
    2042             : 
    2043           0 :       nsize = SIZE(CYCLE)
    2044             : 
    2045             :       ! traverse the path
    2046           0 :       area(:) = 0.0_dp
    2047           0 :       DO ia = 1, nsize
    2048             :          ! contributions from all bead pairs of the current atom
    2049           0 :          DO ib = 1, helium%beads - 1
    2050           0 :             r1(:) = helium%pos(:, CYCLE(ia), ib)
    2051           0 :             r2(:) = helium%pos(:, CYCLE(ia), ib + 1)
    2052           0 :             r12(:) = r2(:) - r1(:)
    2053           0 :             CALL helium_pbc(helium, r1)
    2054           0 :             CALL helium_pbc(helium, r12)
    2055           0 :             r2(:) = r1(:) + r12(:)
    2056           0 :             rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
    2057           0 :             rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
    2058           0 :             rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
    2059           0 :             area(:) = area(:) + rcur(:)
    2060             :          END DO
    2061             :          ! contribution from the last bead of the current atom
    2062             :          ! and the first bead of the next atom
    2063           0 :          i1 = CYCLE(ia)
    2064           0 :          IF (ia .EQ. nsize) THEN
    2065           0 :             i2 = CYCLE(1)
    2066             :          ELSE
    2067           0 :             i2 = CYCLE(ia + 1)
    2068             :          END IF
    2069           0 :          r1(:) = helium%pos(:, i1, helium%beads)
    2070           0 :          r2(:) = helium%pos(:, i2, 1)
    2071           0 :          r12(:) = r2(:) - r1(:)
    2072           0 :          CALL helium_pbc(helium, r1)
    2073           0 :          CALL helium_pbc(helium, r12)
    2074           0 :          r2(:) = r1(:) + r12(:)
    2075           0 :          rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
    2076           0 :          rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
    2077           0 :          rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
    2078           0 :          area(:) = area(:) + rcur(:)
    2079             :       END DO
    2080           0 :       area(:) = 0.5_dp*area(:)
    2081             : 
    2082           0 :    END FUNCTION helium_cycle_projected_area_pbc
    2083             : 
    2084             : ! ***************************************************************************
    2085             : !> \brief  Return total projected area (sum over all cycles)
    2086             : !> \param helium ...
    2087             : !> \return ...
    2088             : !> \date   2014-04-24
    2089             : !> \author Lukasz Walewski
    2090             : !> \note   mostly for sanity checks
    2091             : ! **************************************************************************************************
    2092           0 :    FUNCTION helium_total_projected_area_cyclewise(helium) RESULT(area)
    2093             : 
    2094             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2095             :       REAL(KIND=dp), DIMENSION(3)                        :: area
    2096             : 
    2097             :       INTEGER                                            :: ic
    2098             :       REAL(KIND=dp), DIMENSION(3)                        :: pa
    2099           0 :       TYPE(int_arr_ptr), DIMENSION(:), POINTER           :: cycles
    2100             : 
    2101             : ! decompose the current permutation state into permutation cycles
    2102             : 
    2103           0 :       NULLIFY (cycles)
    2104           0 :       cycles => helium_calc_cycles(helium%permutation)
    2105             : 
    2106           0 :       area(:) = 0.0_dp
    2107           0 :       DO ic = 1, SIZE(cycles)
    2108           0 :          pa = helium_cycle_projected_area(helium, cycles(ic)%iap)
    2109           0 :          area(:) = area(:) + pa(:)
    2110             :       END DO
    2111             : 
    2112           0 :    END FUNCTION helium_total_projected_area_cyclewise
    2113             : 
    2114             : ! ***************************************************************************
    2115             : !> \brief  Return total moment of inertia divided by m_He
    2116             : !> \param helium ...
    2117             : !> \return ...
    2118             : !> \date   2014-04-24
    2119             : !> \author Lukasz Walewski
    2120             : ! **************************************************************************************************
    2121         153 :    FUNCTION helium_total_moment_of_inertia(helium) RESULT(moit)
    2122             : 
    2123             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2124             :       REAL(KIND=dp), DIMENSION(3)                        :: moit
    2125             : 
    2126             :       INTEGER                                            :: ia, ib
    2127             :       REAL(KIND=dp), DIMENSION(3)                        :: com, r1, r12, r2, rcur
    2128             : 
    2129         612 :       com(:) = helium%center(:)
    2130             : 
    2131         612 :       moit(:) = 0.0_dp
    2132        3369 :       DO ia = 1, helium%atoms
    2133             :          ! contributions from all the links of the current atom
    2134       64128 :          DO ib = 1, helium%beads - 1
    2135      243648 :             r1(:) = helium%pos(:, ia, ib) - com(:)
    2136      243648 :             r2(:) = helium%pos(:, ia, ib + 1) - com(:)
    2137             :             ! comment out for non-PBC version -->
    2138      243648 :             r12(:) = r2(:) - r1(:)
    2139       60912 :             CALL helium_pbc(helium, r1)
    2140       60912 :             CALL helium_pbc(helium, r12)
    2141      243648 :             r2(:) = r1(:) + r12(:)
    2142             :             ! comment out for non-PBC version <--
    2143       60912 :             rcur(1) = r1(2)*r2(2) + r1(3)*r2(3)
    2144       60912 :             rcur(2) = r1(3)*r2(3) + r1(1)*r2(1)
    2145       60912 :             rcur(3) = r1(1)*r2(1) + r1(2)*r2(2)
    2146      246864 :             moit(:) = moit(:) + rcur(:)
    2147             :          END DO
    2148             :          ! contribution from the last bead of the current atom
    2149             :          ! and the first bead of the next atom
    2150       12864 :          r1(:) = helium%pos(:, ia, helium%beads) - com(:)
    2151       12864 :          r2(:) = helium%pos(:, helium%permutation(ia), 1) - com(:)
    2152             :          ! comment out for non-PBC version -->
    2153       12864 :          r12(:) = r2(:) - r1(:)
    2154        3216 :          CALL helium_pbc(helium, r1)
    2155        3216 :          CALL helium_pbc(helium, r12)
    2156       12864 :          r2(:) = r1(:) + r12(:)
    2157             :          ! comment out for non-PBC version <--
    2158        3216 :          rcur(1) = r1(2)*r2(2) + r1(3)*r2(3)
    2159        3216 :          rcur(2) = r1(3)*r2(3) + r1(1)*r2(1)
    2160        3216 :          rcur(3) = r1(1)*r2(1) + r1(2)*r2(2)
    2161       13017 :          moit(:) = moit(:) + rcur(:)
    2162             :       END DO
    2163         612 :       moit(:) = moit(:)/helium%beads
    2164             : 
    2165         153 :    END FUNCTION helium_total_moment_of_inertia
    2166             : 
    2167             : ! ***************************************************************************
    2168             : !> \brief  Return link moment of inertia divided by m_He
    2169             : !> \param helium ...
    2170             : !> \param ia ...
    2171             : !> \param ib ...
    2172             : !> \return ...
    2173             : !> \date   2014-04-24
    2174             : !> \author Lukasz Walewski
    2175             : ! **************************************************************************************************
    2176           0 :    FUNCTION helium_link_moment_of_inertia(helium, ia, ib) RESULT(moit)
    2177             : 
    2178             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2179             :       INTEGER, INTENT(IN)                                :: ia, ib
    2180             :       REAL(KIND=dp), DIMENSION(3)                        :: moit
    2181             : 
    2182             :       INTEGER                                            :: ja1, ja2, jb1, jb2
    2183             :       REAL(KIND=dp), DIMENSION(3)                        :: com, r1, r12, r2
    2184             : 
    2185           0 :       com(:) = helium%center(:)
    2186             : 
    2187           0 :       IF (ib .EQ. helium%beads) THEN
    2188           0 :          ja1 = ia
    2189           0 :          ja2 = helium%permutation(ia)
    2190           0 :          jb1 = ib
    2191           0 :          jb2 = 1
    2192             :       ELSE
    2193           0 :          ja1 = ia
    2194           0 :          ja2 = ia
    2195           0 :          jb1 = ib
    2196           0 :          jb2 = ib + 1
    2197             :       END IF
    2198           0 :       r1(:) = helium%pos(:, ja1, jb1) - com(:)
    2199           0 :       r2(:) = helium%pos(:, ja2, jb2) - com(:)
    2200             :       ! comment out for non-PBC version -->
    2201           0 :       r12(:) = r2(:) - r1(:)
    2202           0 :       CALL helium_pbc(helium, r1)
    2203           0 :       CALL helium_pbc(helium, r12)
    2204           0 :       r2(:) = r1(:) + r12(:)
    2205             :       ! comment out for non-PBC version <--
    2206           0 :       moit(1) = r1(2)*r2(2) + r1(3)*r2(3)
    2207           0 :       moit(2) = r1(3)*r2(3) + r1(1)*r2(1)
    2208           0 :       moit(3) = r1(1)*r2(1) + r1(2)*r2(2)
    2209           0 :       moit(:) = moit(:)/helium%beads
    2210             : 
    2211           0 :    END FUNCTION helium_link_moment_of_inertia
    2212             : 
    2213             : ! ***************************************************************************
    2214             : !> \brief  Return total moment of inertia (sum over all links)
    2215             : !> \param helium ...
    2216             : !> \return ...
    2217             : !> \date   2014-04-24
    2218             : !> \author Lukasz Walewski
    2219             : !> \note   mostly for sanity checks
    2220             : ! **************************************************************************************************
    2221           0 :    FUNCTION helium_total_moment_of_inertia_linkwise(helium) RESULT(moit)
    2222             : 
    2223             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2224             :       REAL(KIND=dp), DIMENSION(3)                        :: moit
    2225             : 
    2226             :       INTEGER                                            :: ia, ib
    2227             : 
    2228           0 :       moit(:) = 0.0_dp
    2229           0 :       DO ia = 1, helium%atoms
    2230           0 :          DO ib = 1, helium%beads
    2231           0 :             moit(:) = moit(:) + helium_link_moment_of_inertia(helium, ia, ib)
    2232             :          END DO
    2233             :       END DO
    2234             : 
    2235           0 :    END FUNCTION helium_total_moment_of_inertia_linkwise
    2236             : 
    2237             : ! ***************************************************************************
    2238             : !> \brief  Return moment of inertia of a cycle divided by m_He
    2239             : !> \param helium ...
    2240             : !> \param CYCLE ...
    2241             : !> \param pos ...
    2242             : !> \return ...
    2243             : !> \date   2014-04-24
    2244             : !> \author Lukasz Walewski
    2245             : ! **************************************************************************************************
    2246           0 :    FUNCTION helium_cycle_moment_of_inertia(helium, CYCLE, pos) RESULT(moit)
    2247             : 
    2248             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2249             :       INTEGER, DIMENSION(:), POINTER                     :: CYCLE
    2250             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos
    2251             :       REAL(KIND=dp), DIMENSION(3)                        :: moit
    2252             : 
    2253             :       INTEGER                                            :: i1, i2, ia, ib, nsize
    2254             :       REAL(KIND=dp), DIMENSION(3)                        :: com, rcur, ri, rj
    2255             : 
    2256           0 :       nsize = SIZE(CYCLE)
    2257             : 
    2258             :       ! traverse the path
    2259           0 :       moit(:) = 0.0_dp
    2260           0 :       com(:) = helium_com(helium)
    2261           0 :       DO ia = 1, nsize
    2262             :          ! contributions from all bead pairs of the current atom
    2263           0 :          DO ib = 1, helium%beads - 1
    2264           0 :             ri = pos(:, CYCLE(ia), ib) - com(:)
    2265           0 :             rj = pos(:, CYCLE(ia), ib + 1) - com(:)
    2266           0 :             rcur(1) = ri(2)*rj(2) + ri(3)*rj(3)
    2267           0 :             rcur(2) = ri(3)*rj(3) + ri(1)*rj(1)
    2268           0 :             rcur(3) = ri(1)*rj(1) + ri(2)*rj(2)
    2269           0 :             moit(:) = moit(:) + rcur(:)
    2270             :          END DO
    2271             :          ! contribution from the last bead of the current atom
    2272             :          ! and the first bead of the next atom
    2273           0 :          i1 = CYCLE(ia)
    2274           0 :          IF (ia .EQ. nsize) THEN
    2275           0 :             i2 = CYCLE(1)
    2276             :          ELSE
    2277           0 :             i2 = CYCLE(ia + 1)
    2278             :          END IF
    2279             :          ! rotation invariant bead index
    2280           0 :          ri = pos(:, i1, helium%beads) - com(:)
    2281           0 :          rj = pos(:, i2, 1) - com(:)
    2282           0 :          rcur(1) = ri(2)*rj(2) + ri(3)*rj(3)
    2283           0 :          rcur(2) = ri(3)*rj(3) + ri(1)*rj(1)
    2284           0 :          rcur(3) = ri(1)*rj(1) + ri(2)*rj(2)
    2285           0 :          moit(:) = moit(:) + rcur(:)
    2286             :       END DO
    2287           0 :       moit(:) = moit(:)/helium%beads
    2288             : 
    2289           0 :    END FUNCTION helium_cycle_moment_of_inertia
    2290             : 
    2291             : ! ***************************************************************************
    2292             : !> \brief  Return total moment of inertia (sum over all cycles)
    2293             : !> \param helium ...
    2294             : !> \return ...
    2295             : !> \date   2014-04-24
    2296             : !> \author Lukasz Walewski
    2297             : !> \note   mostly for sanity checks
    2298             : ! **************************************************************************************************
    2299           0 :    FUNCTION helium_total_moment_of_inertia_cyclewise(helium) RESULT(moit)
    2300             : 
    2301             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2302             :       REAL(KIND=dp), DIMENSION(3)                        :: moit
    2303             : 
    2304             :       INTEGER                                            :: ic
    2305             :       REAL(KIND=dp), DIMENSION(3)                        :: pa
    2306           0 :       TYPE(int_arr_ptr), DIMENSION(:), POINTER           :: cycles
    2307             : 
    2308             : ! decompose the current permutation state into permutation cycles
    2309             : 
    2310           0 :       NULLIFY (cycles)
    2311           0 :       cycles => helium_calc_cycles(helium%permutation)
    2312             : 
    2313           0 :       moit(:) = 0.0_dp
    2314           0 :       DO ic = 1, SIZE(cycles)
    2315           0 :          pa = helium_cycle_moment_of_inertia(helium, cycles(ic)%iap, helium%pos)
    2316           0 :          moit(:) = moit(:) + pa(:)
    2317             :       END DO
    2318             : 
    2319           0 :       DEALLOCATE (cycles)
    2320             : 
    2321           0 :    END FUNCTION helium_total_moment_of_inertia_cyclewise
    2322             : 
    2323             : ! ***************************************************************************
    2324             : !> \brief  Set coordinate system, e.g. for RHO calculations
    2325             : !> \param helium ...
    2326             : !> \param pint_env ...
    2327             : !> \date   2014-04-25
    2328             : !> \author Lukasz Walewski
    2329             : !> \note   Sets the origin (center of the coordinate system) wrt which
    2330             : !>         spatial distribution functions are calculated.
    2331             : !> \note   It can be extended later to set the axes of the coordinate system
    2332             : !>         as well, e.g. for dynamic analysis with moving solute.
    2333             : ! **************************************************************************************************
    2334          27 :    PURE SUBROUTINE helium_update_coord_system(helium, pint_env)
    2335             : 
    2336             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    2337             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2338             : 
    2339          27 :       IF (helium%solute_present) THEN
    2340          68 :          helium%center(:) = pint_com_pos(pint_env)
    2341             :       ELSE
    2342          10 :          IF (helium%periodic) THEN
    2343          40 :             helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
    2344             :          ELSE
    2345           0 :             helium%center(:) = helium_com(helium)
    2346             :          END IF
    2347             :       END IF
    2348             : 
    2349          27 :    END SUBROUTINE helium_update_coord_system
    2350             : 
    2351             : ! ***************************************************************************
    2352             : !> \brief  Set coordinate system for RDF calculations
    2353             : !> \param helium ...
    2354             : !> \param pint_env ...
    2355             : !> \date   2014-04-25
    2356             : !> \par    History
    2357             : !>         2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
    2358             : !> \author Lukasz Walewski
    2359             : !> \note   Sets the centers wrt which radial distribution functions are
    2360             : !>         calculated.
    2361             : ! **************************************************************************************************
    2362          12 :    PURE SUBROUTINE helium_set_rdf_coord_system(helium, pint_env)
    2363             : 
    2364             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    2365             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2366             : 
    2367             :       INTEGER                                            :: i, j
    2368             : 
    2369          12 :       IF (helium%solute_present .AND. helium%rdf_sol_he) THEN
    2370             :          ! Account for unequal number of beads for solute and helium
    2371           0 :          DO i = 1, helium%beads
    2372           0 :             j = ((i - 1)*helium%solute_beads)/helium%beads + 1
    2373           0 :             helium%rdf_centers(i, :) = pint_env%x(j, :)
    2374             :          END DO
    2375             :       END IF
    2376             : 
    2377          12 :    END SUBROUTINE helium_set_rdf_coord_system
    2378             : 
    2379             : ! ***************************************************************************
    2380             : !> \brief  Calculate center of mass
    2381             : !> \param helium ...
    2382             : !> \return ...
    2383             : !> \date   2014-04-09
    2384             : !> \author Lukasz Walewski
    2385             : ! **************************************************************************************************
    2386           0 :    PURE FUNCTION helium_com(helium) RESULT(com)
    2387             : 
    2388             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2389             :       REAL(KIND=dp), DIMENSION(3)                        :: com
    2390             : 
    2391             :       INTEGER                                            :: ia, ib
    2392             : 
    2393           0 :       com(:) = 0.0_dp
    2394           0 :       DO ia = 1, helium%atoms
    2395           0 :          DO ib = 1, helium%beads
    2396           0 :             com(:) = com(:) + helium%pos(:, ia, ib)
    2397             :          END DO
    2398             :       END DO
    2399           0 :       com(:) = com(:)/helium%atoms/helium%beads
    2400             : 
    2401           0 :    END FUNCTION helium_com
    2402             : 
    2403             : ! ***************************************************************************
    2404             : !> \brief  Return link vector, i.e. displacement vector of two consecutive
    2405             : !>         beads along the path starting at bead ib of atom ia
    2406             : !> \param helium ...
    2407             : !> \param ia ...
    2408             : !> \param ib ...
    2409             : !> \return ...
    2410             : !> \author Lukasz Walewski
    2411             : ! **************************************************************************************************
    2412           0 :    FUNCTION helium_link_vector(helium, ia, ib) RESULT(lvec)
    2413             : 
    2414             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2415             :       INTEGER, INTENT(IN)                                :: ia, ib
    2416             :       REAL(KIND=dp), DIMENSION(3)                        :: lvec
    2417             : 
    2418             :       INTEGER                                            :: ia1, ia2, ib1, ib2
    2419           0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: r1, r2
    2420             : 
    2421           0 :       IF (ib .EQ. helium%beads) THEN
    2422           0 :          ia1 = ia
    2423           0 :          ia2 = helium%permutation(ia)
    2424           0 :          ib1 = ib
    2425           0 :          ib2 = 1
    2426             :       ELSE
    2427           0 :          ia1 = ia
    2428           0 :          ia2 = ia
    2429           0 :          ib1 = ib
    2430           0 :          ib2 = ib + 1
    2431             :       END IF
    2432           0 :       r1 => helium%pos(:, ia1, ib1)
    2433           0 :       r2 => helium%pos(:, ia2, ib2)
    2434           0 :       lvec(:) = r2(:) - r1(:)
    2435           0 :       CALL helium_pbc(helium, lvec)
    2436             : 
    2437           0 :    END FUNCTION helium_link_vector
    2438             : 
    2439             : ! **************************************************************************************************
    2440             : !> \brief ...
    2441             : !> \param helium ...
    2442             : !> \param ia ...
    2443             : !> \param ib ...
    2444             : !> \return ...
    2445             : ! **************************************************************************************************
    2446           0 :    PURE FUNCTION helium_rperpendicular(helium, ia, ib) RESULT(rperp)
    2447             : 
    2448             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2449             :       INTEGER, INTENT(IN)                                :: ia, ib
    2450             :       REAL(KIND=dp), DIMENSION(3)                        :: rperp
    2451             : 
    2452             :       ASSOCIATE (ri => helium%pos(:, ia, ib))
    2453           0 :          rperp(1) = SQRT(ri(2)*ri(2) + ri(3)*ri(3))
    2454           0 :          rperp(2) = SQRT(ri(3)*ri(3) + ri(1)*ri(1))
    2455           0 :          rperp(3) = SQRT(ri(1)*ri(1) + ri(2)*ri(2))
    2456             :       END ASSOCIATE
    2457             : 
    2458           0 :    END FUNCTION helium_rperpendicular
    2459             : 
    2460             : ! ***************************************************************************
    2461             : !> \brief  Convert a winding number vector from real number representation
    2462             : !>         (in internal units) to integer number representation (in cell
    2463             : !>         vector units)
    2464             : !> \param helium ...
    2465             : !> \param wnum ...
    2466             : !> \return ...
    2467             : !> \author Lukasz Walewski
    2468             : ! **************************************************************************************************
    2469         250 :    FUNCTION helium_wnumber_to_integer(helium, wnum) RESULT(inum)
    2470             : 
    2471             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2472             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: wnum
    2473             :       INTEGER, DIMENSION(3)                              :: inum
    2474             : 
    2475             :       REAL(KIND=dp), DIMENSION(3)                        :: wcur
    2476             : 
    2477         250 :       CALL DGEMV('N', 3, 3, 1.0_dp, helium%cell_m_inv, SIZE(helium%cell_m_inv, 1), wnum, 1, 0.0_dp, wcur, 1)
    2478        1000 :       inum(:) = NINT(wcur(:))
    2479             : 
    2480         250 :    END FUNCTION helium_wnumber_to_integer
    2481             : 
    2482             : ! ***************************************************************************
    2483             : !> \brief  Given the atom index and permutation state returns .TRUE. if the
    2484             : !>         atom belongs to a winding path, .FASLE. otherwise.
    2485             : !> \param helium ...
    2486             : !> \param atmidx ...
    2487             : !> \param pos ...
    2488             : !> \param permutation ...
    2489             : !> \return ...
    2490             : !> \date   2010-09-21
    2491             : !> \author Lukasz Walewski
    2492             : !> \note   The path winds around the periodic box if any component of its
    2493             : !>         widing number vector differs from 0.
    2494             : ! **************************************************************************************************
    2495         250 :    FUNCTION helium_is_winding(helium, atmidx, pos, permutation) RESULT(is_winding)
    2496             : 
    2497             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
    2498             :       INTEGER, INTENT(IN)                                :: atmidx
    2499             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos
    2500             :       INTEGER, DIMENSION(:), POINTER                     :: permutation
    2501             :       LOGICAL                                            :: is_winding
    2502             : 
    2503             :       INTEGER                                            :: ic
    2504             :       INTEGER, DIMENSION(3)                              :: inum
    2505             :       INTEGER, DIMENSION(:), POINTER                     :: CYCLE
    2506             :       REAL(KIND=dp), DIMENSION(3)                        :: wnum
    2507             : 
    2508         250 :       is_winding = .FALSE.
    2509             :       NULLIFY (CYCLE)
    2510         250 :       CYCLE => helium_cycle_of(atmidx, permutation)
    2511        1250 :       wnum(:) = bohr*helium_cycle_winding_number(helium, CYCLE, pos)
    2512         250 :       inum(:) = helium_wnumber_to_integer(helium, wnum)
    2513        1000 :       DO ic = 1, 3
    2514        1000 :          IF (ABS(inum(ic)) .GT. 0) THEN
    2515             :             is_winding = .TRUE.
    2516             :             EXIT
    2517             :          END IF
    2518             :       END DO
    2519         250 :       DEALLOCATE (CYCLE)
    2520             : 
    2521         250 :    END FUNCTION helium_is_winding
    2522             : 
    2523             : END MODULE helium_common

Generated by: LCOV version 1.15