LCOV - code coverage report
Current view: top level - src/motion - helium_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 1081 1330 81.3 %
Date: 2023-03-30 11:55:16 Functions: 22 22 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief  Methods dealing with helium_solvent_type
      10             : !> \author Lukasz Walewski
      11             : !> \date   2009-06-10
      12             : ! **************************************************************************************************
      13             : MODULE helium_methods
      14             : 
      15             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      16             :    USE bibliography,                    ONLY: Walewski2014,&
      17             :                                               cite_reference
      18             :    USE cell_types,                      ONLY: cell_type,&
      19             :                                               get_cell
      20             :    USE cp_files,                        ONLY: close_file,&
      21             :                                               open_file
      22             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23             :                                               cp_logger_type
      24             :    USE cp_output_handling,              ONLY: cp_printkey_is_on
      25             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      26             :                                               cp_subsys_type
      27             :    USE f77_interface,                   ONLY: f_env_add_defaults,&
      28             :                                               f_env_rm_defaults,&
      29             :                                               f_env_type
      30             :    USE force_env_types,                 ONLY: force_env_get
      31             :    USE helium_common,                   ONLY: helium_com,&
      32             :                                               helium_pbc
      33             :    USE helium_interactions,             ONLY: helium_vij
      34             :    USE helium_io,                       ONLY: helium_write_line,&
      35             :                                               helium_write_setup
      36             :    USE helium_sampling,                 ONLY: helium_sample
      37             :    USE helium_types,                    ONLY: helium_solvent_p_type,&
      38             :                                               helium_solvent_type,&
      39             :                                               rho_atom_number,&
      40             :                                               rho_moment_of_inertia,&
      41             :                                               rho_num,&
      42             :                                               rho_projected_area,&
      43             :                                               rho_winding_cycle,&
      44             :                                               rho_winding_number
      45             :    USE input_constants,                 ONLY: helium_cell_shape_cube,&
      46             :                                               helium_cell_shape_octahedron,&
      47             :                                               helium_sampling_ceperley,&
      48             :                                               helium_sampling_worm,&
      49             :                                               helium_solute_intpot_none
      50             :    USE input_section_types,             ONLY: section_vals_get,&
      51             :                                               section_vals_get_subs_vals,&
      52             :                                               section_vals_type,&
      53             :                                               section_vals_val_get,&
      54             :                                               section_vals_val_set
      55             :    USE kinds,                           ONLY: default_path_length,&
      56             :                                               default_string_length,&
      57             :                                               dp,&
      58             :                                               max_line_length
      59             :    USE mathconstants,                   ONLY: pi,&
      60             :                                               twopi
      61             :    USE message_passing,                 ONLY: mp_comm_type
      62             :    USE parallel_rng_types,              ONLY: GAUSSIAN,&
      63             :                                               UNIFORM,&
      64             :                                               rng_stream_p_type,&
      65             :                                               rng_stream_type
      66             :    USE particle_list_types,             ONLY: particle_list_type
      67             :    USE physcon,                         ONLY: a_mass,&
      68             :                                               angstrom,&
      69             :                                               boltzmann,&
      70             :                                               h_bar,&
      71             :                                               kelvin,&
      72             :                                               massunit
      73             :    USE pint_public,                     ONLY: pint_com_pos
      74             :    USE pint_types,                      ONLY: pint_env_type
      75             :    USE splines_methods,                 ONLY: init_spline,&
      76             :                                               init_splinexy
      77             :    USE splines_types,                   ONLY: spline_data_create,&
      78             :                                               spline_data_release
      79             : #include "../base/base_uses.f90"
      80             : 
      81             :    IMPLICIT NONE
      82             : 
      83             :    PRIVATE
      84             : 
      85             :    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
      86             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_methods'
      87             : 
      88             :    PUBLIC :: helium_create
      89             :    PUBLIC :: helium_init
      90             :    PUBLIC :: helium_release
      91             : 
      92             : CONTAINS
      93             : 
      94             : ! ***************************************************************************
      95             : !> \brief  Data-structure that holds all needed information about
      96             : !>         (superfluid) helium solvent
      97             : !> \param helium_env ...
      98             : !> \param input ...
      99             : !> \param solute ...
     100             : !> \par    History
     101             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     102             : !> \author hforbert
     103             : ! **************************************************************************************************
     104          24 :    SUBROUTINE helium_create(helium_env, input, solute)
     105             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     106             :       TYPE(section_vals_type), POINTER                   :: input
     107             :       TYPE(pint_env_type), INTENT(IN), OPTIONAL          :: solute
     108             : 
     109             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_create'
     110             : 
     111             :       CHARACTER(len=default_path_length)                 :: msg_str, potential_file_name
     112             :       INTEGER                                            :: color_sub, handle, i, input_unit, isize, &
     113             :                                                             itmp, j, k, mepos, nlines, ntab, &
     114             :                                                             num_env, pdx
     115          24 :       INTEGER, DIMENSION(:), POINTER                     :: env_all
     116             :       LOGICAL                                            :: expl_cell, expl_dens, expl_nats, &
     117             :                                                             expl_pot, explicit, ltmp
     118             :       REAL(KIND=dp)                                      :: cgeof, dx, he_mass, mHe, rtmp, T, tau, &
     119             :                                                             tcheck, x1, x_spline
     120          24 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: pot_transfer
     121             :       TYPE(cp_logger_type), POINTER                      :: logger
     122             :       TYPE(mp_comm_type)                                 :: new_comm
     123             :       TYPE(section_vals_type), POINTER                   :: helium_section, input_worm
     124             : 
     125          24 :       CALL timeset(routineN, handle)
     126             : 
     127          24 :       CALL cite_reference(Walewski2014)
     128          24 :       NULLIFY (logger)
     129          24 :       logger => cp_get_default_logger()
     130             : 
     131          24 :       NULLIFY (helium_section)
     132             :       helium_section => section_vals_get_subs_vals(input, &
     133          24 :                                                    "MOTION%PINT%HELIUM")
     134          24 :       CALL section_vals_get(helium_section, explicit=explicit)
     135          24 :       CPASSERT(explicit)
     136             : 
     137             :       ! get number of environments
     138             :       CALL section_vals_val_get(helium_section, "NUM_ENV", &
     139          24 :                                 explicit=explicit)
     140          24 :       IF (explicit) THEN
     141             :          CALL section_vals_val_get(helium_section, "NUM_ENV", &
     142          24 :                                    i_val=num_env)
     143             :       ELSE
     144           0 :          num_env = logger%para_env%num_pe
     145             :       END IF
     146          24 :       CPASSERT(num_env >= 0)
     147          24 :       IF (num_env .NE. logger%para_env%num_pe) THEN
     148           4 :          msg_str = "NUM_ENV not equal to number of processors"
     149           4 :          CPWARN(msg_str)
     150             :       END IF
     151             :       ! calculate number of tasks for each processor
     152             :       mepos = num_env/logger%para_env%num_pe &
     153          24 :               + MIN(MOD(num_env, logger%para_env%num_pe)/(logger%para_env%mepos + 1), 1)
     154             :       ! gather result
     155          24 :       NULLIFY (env_all)
     156          72 :       ALLOCATE (env_all(logger%para_env%num_pe))
     157          72 :       env_all(:) = 0
     158          72 :       CALL logger%para_env%allgather(mepos, env_all)
     159             : 
     160             :       ! create new communicator for processors with helium_env
     161          24 :       IF (mepos == 0) THEN
     162           0 :          color_sub = 0
     163             :       ELSE
     164          24 :          color_sub = 1
     165             :       END IF
     166          24 :       CALL new_comm%from_split(logger%para_env, color_sub)
     167             :       ! release new_comm for processors without helium_env
     168          24 :       IF (mepos == 0) THEN
     169           0 :          CALL new_comm%free()
     170             :       END IF
     171             : 
     172          24 :       NULLIFY (helium_env)
     173          24 :       IF (mepos .GT. 0) THEN
     174         102 :          ALLOCATE (helium_env(mepos))
     175          54 :          DO k = 1, mepos
     176          30 :             helium_env(k)%comm = new_comm
     177          30 :             NULLIFY (helium_env(k)%env_all)
     178          30 :             helium_env(k)%env_all => env_all
     179          30 :             ALLOCATE (helium_env(k)%helium)
     180          30 :             NULLIFY (helium_env(k)%helium%input)
     181          30 :             helium_env(k)%helium%input => input
     182          54 :             helium_env(k)%helium%num_env = num_env
     183             :          END DO
     184             :          ! RNG state create & init
     185          24 :          CALL helium_rng_init(helium_env)
     186             : 
     187          54 :          DO k = 1, mepos
     188             :             NULLIFY (helium_env(k)%helium%ptable, &
     189          30 :                      helium_env(k)%helium%permutation, &
     190          30 :                      helium_env(k)%helium%savepermutation, &
     191          30 :                      helium_env(k)%helium%iperm, &
     192          30 :                      helium_env(k)%helium%saveiperm, &
     193          30 :                      helium_env(k)%helium%itmp_atoms_1d, &
     194          30 :                      helium_env(k)%helium%ltmp_atoms_1d, &
     195          30 :                      helium_env(k)%helium%itmp_atoms_np_1d, &
     196          30 :                      helium_env(k)%helium%pos, &
     197          30 :                      helium_env(k)%helium%savepos, &
     198          30 :                      helium_env(k)%helium%work, &
     199          30 :                      helium_env(k)%helium%force_avrg, &
     200          30 :                      helium_env(k)%helium%force_inst, &
     201          30 :                      helium_env(k)%helium%rtmp_3_np_1d, &
     202          30 :                      helium_env(k)%helium%rtmp_p_ndim_1d, &
     203          30 :                      helium_env(k)%helium%rtmp_p_ndim_np_1d, &
     204          30 :                      helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
     205          30 :                      helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
     206          30 :                      helium_env(k)%helium%rtmp_p_ndim_2d, &
     207          30 :                      helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
     208          30 :                      helium_env(k)%helium%tmatrix, helium_env(k)%helium%pmatrix, &
     209          30 :                      helium_env(k)%helium%nmatrix, helium_env(k)%helium%ipmatrix, &
     210          30 :                      helium_env(k)%helium%u0, helium_env(k)%helium%e0, &
     211          30 :                      helium_env(k)%helium%uoffdiag, helium_env(k)%helium%eoffdiag, &
     212          30 :                      helium_env(k)%helium%vij, &
     213          30 :                      helium_env(k)%helium%rdf_inst, &
     214          30 :                      helium_env(k)%helium%plength_avrg, &
     215          30 :                      helium_env(k)%helium%plength_inst, &
     216          30 :                      helium_env(k)%helium%atom_plength, &
     217          30 :                      helium_env(k)%helium%ename &
     218          30 :                      )
     219             : 
     220          30 :             helium_env(k)%helium%accepts = 0
     221          30 :             helium_env(k)%helium%relrot = 0
     222             : 
     223             :             ! check if solute is present in our simulation
     224          30 :             helium_env(k)%helium%solute_present = .FALSE.
     225          30 :             helium_env(k)%helium%solute_atoms = 0
     226          30 :             helium_env(k)%helium%solute_beads = 0
     227             :             CALL section_vals_val_get( &
     228             :                helium_section, &
     229             :                "HELIUM_ONLY", &
     230          30 :                l_val=ltmp)
     231          30 :             IF (.NOT. ltmp) THEN
     232          20 :                IF (PRESENT(solute)) THEN
     233          20 :                   helium_env(k)%helium%solute_present = .TRUE.
     234          20 :                   helium_env(k)%helium%solute_atoms = solute%ndim/3
     235          20 :                   helium_env(k)%helium%solute_beads = solute%p
     236             :                END IF
     237             :             END IF
     238             : 
     239             :             CALL section_vals_val_get(helium_section, "NBEADS", &
     240          30 :                                       i_val=helium_env(k)%helium%beads)
     241             :             CALL section_vals_val_get(helium_section, "INOROT", &
     242          30 :                                       i_val=helium_env(k)%helium%iter_norot)
     243             :             CALL section_vals_val_get(helium_section, "IROT", &
     244          30 :                                       i_val=helium_env(k)%helium%iter_rot)
     245             : 
     246             :             ! get number of steps and current step number from PINT
     247             :             CALL section_vals_val_get(input, "MOTION%PINT%ITERATION", &
     248          30 :                                       i_val=itmp)
     249          30 :             helium_env(k)%helium%first_step = itmp
     250             :             CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
     251          30 :                                       explicit=explicit)
     252          30 :             IF (explicit) THEN
     253             :                CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
     254           0 :                                          i_val=itmp)
     255           0 :                helium_env(k)%helium%last_step = itmp
     256             :                helium_env(k)%helium%num_steps = helium_env(k)%helium%last_step &
     257           0 :                                                 - helium_env(k)%helium%first_step
     258             :             ELSE
     259             :                CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", &
     260          30 :                                          i_val=itmp)
     261          30 :                helium_env(k)%helium%num_steps = itmp
     262             :                helium_env(k)%helium%last_step = helium_env(k)%helium%first_step &
     263          30 :                                                 + helium_env(k)%helium%num_steps
     264             :             END IF
     265             : 
     266             :             ! boundary conditions
     267             :             CALL section_vals_val_get(helium_section, "PERIODIC", &
     268          30 :                                       l_val=helium_env(k)%helium%periodic)
     269             :             CALL section_vals_val_get(helium_section, "CELL_SHAPE", &
     270          30 :                                       i_val=helium_env(k)%helium%cell_shape)
     271             : 
     272             :             CALL section_vals_val_get(helium_section, "DROPLET_RADIUS", &
     273          30 :                                       r_val=helium_env(k)%helium%droplet_radius)
     274             : 
     275             :             ! Set density Rho, number of atoms N and volume V ( Rho = N / V ).
     276             :             ! Allow only 2 out of 3 values to be defined at the same time, calculate
     277             :             ! the third.
     278             :             ! Note, that DENSITY and NATOMS keywords have default values, while
     279             :             ! CELL_SIZE does not. Thus if CELL_SIZE is given explicitly then one and
     280             :             ! only one of the two remaining options must be give explicitly as well.
     281             :             ! If CELL_SIZE is not given explicitly then all four combinations of the
     282             :             ! two other options are valid.
     283             :             CALL section_vals_val_get(helium_section, "DENSITY", &
     284          30 :                                       explicit=expl_dens, r_val=helium_env(k)%helium%density)
     285             :             CALL section_vals_val_get(helium_section, "NATOMS", &
     286          30 :                                       explicit=expl_nats, i_val=helium_env(k)%helium%atoms)
     287             :             CALL section_vals_val_get(helium_section, "CELL_SIZE", &
     288          30 :                                       explicit=expl_cell)
     289          30 :             cgeof = 1.0_dp
     290          30 :             IF (helium_env(k)%helium%periodic) THEN
     291          16 :                IF (helium_env(k)%helium%cell_shape .EQ. helium_cell_shape_octahedron) cgeof = 2.0_dp
     292             :             END IF
     293          30 :             rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density)**(1.0_dp/3.0_dp)
     294          30 :             IF (.NOT. expl_cell) THEN
     295          30 :                helium_env(k)%helium%cell_size = rtmp
     296             :             ELSE
     297             :                CALL section_vals_val_get(helium_section, "CELL_SIZE", &
     298           0 :                                          r_val=helium_env(k)%helium%cell_size)
     299             :                ! only more work if not all three values are consistent:
     300           0 :                IF (ABS(helium_env(k)%helium%cell_size - rtmp) .GT. 100.0_dp*EPSILON(0.0_dp)* &
     301             :                    (ABS(helium_env(k)%helium%cell_size) + rtmp)) THEN
     302           0 :                   IF (expl_dens .AND. expl_nats) THEN
     303             :                      msg_str = "DENSITY, NATOMS and CELL_SIZE options "// &
     304           0 :                                "contradict each other"
     305           0 :                      CPWARN(msg_str)
     306             :                   END IF
     307             :                   !ok we have enough freedom to resolve the conflict:
     308           0 :                   IF (.NOT. expl_dens) THEN
     309           0 :                      helium_env(k)%helium%density = cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%cell_size**3.0_dp
     310           0 :                      IF (.NOT. expl_nats) THEN
     311             :                         msg_str = "CELL_SIZE defined but neither "// &
     312           0 :                                   "NATOMS nor DENSITY given, using default NATOMS."
     313           0 :                         CPWARN(msg_str)
     314             :                      END IF
     315             :                   ELSE ! ( expl_dens .AND. .NOT. expl_nats )
     316             :                      ! calculate the nearest number of atoms for given conditions
     317             :                      helium_env(k)%helium%atoms = NINT(helium_env(k)%helium%density* &
     318           0 :                                                        helium_env(k)%helium%cell_size**3.0_dp/cgeof)
     319             :                      ! adjust cell size to maintain correct density
     320             :                      ! (should be a small correction)
     321             :                      rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density &
     322           0 :                              )**(1.0_dp/3.0_dp)
     323           0 :                      IF (ABS(helium_env(k)%helium%cell_size - rtmp) .GT. 100.0_dp*EPSILON(0.0_dp) &
     324             :                          *(ABS(helium_env(k)%helium%cell_size) + rtmp)) THEN
     325             :                         msg_str = "Adjusting actual cell size "// &
     326           0 :                                   "to maintain correct density."
     327           0 :                         CPWARN(msg_str)
     328           0 :                         helium_env(k)%helium%cell_size = rtmp
     329             :                      END IF
     330             :                   END IF
     331             :                END IF
     332             :             END IF
     333          30 :             helium_env(k)%helium%cell_size_inv = 1.0_dp/helium_env(k)%helium%cell_size
     334             :             ! From now on helium%density, helium%atoms and helium%cell_size are
     335             :             ! correctly defined.
     336             : 
     337             :             ! set the M matrix for winding number calculations
     338         114 :             SELECT CASE (helium_env(k)%helium%cell_shape)
     339             : 
     340             :             CASE (helium_cell_shape_octahedron)
     341          20 :                helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size
     342          20 :                helium_env(k)%helium%cell_m(2, 1) = 0.0_dp
     343          20 :                helium_env(k)%helium%cell_m(3, 1) = 0.0_dp
     344          20 :                helium_env(k)%helium%cell_m(1, 2) = 0.0_dp
     345          20 :                helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size
     346          20 :                helium_env(k)%helium%cell_m(3, 2) = 0.0_dp
     347          20 :                helium_env(k)%helium%cell_m(1, 3) = helium_env(k)%helium%cell_size/2.0_dp
     348          20 :                helium_env(k)%helium%cell_m(2, 3) = helium_env(k)%helium%cell_size/2.0_dp
     349          20 :                helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size/2.0_dp
     350          20 :                helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size
     351          20 :                helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp
     352          20 :                helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp
     353          20 :                helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp
     354          20 :                helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size
     355          20 :                helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp
     356          20 :                helium_env(k)%helium%cell_m_inv(1, 3) = -1.0_dp/helium_env(k)%helium%cell_size
     357          20 :                helium_env(k)%helium%cell_m_inv(2, 3) = -1.0_dp/helium_env(k)%helium%cell_size
     358          20 :                helium_env(k)%helium%cell_m_inv(3, 3) = 2.0_dp/helium_env(k)%helium%cell_size
     359             :             CASE (helium_cell_shape_cube)
     360             : 
     361          10 :                helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size
     362          10 :                helium_env(k)%helium%cell_m(2, 1) = 0.0_dp
     363          10 :                helium_env(k)%helium%cell_m(3, 1) = 0.0_dp
     364          10 :                helium_env(k)%helium%cell_m(1, 2) = 0.0_dp
     365          10 :                helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size
     366          10 :                helium_env(k)%helium%cell_m(3, 2) = 0.0_dp
     367          10 :                helium_env(k)%helium%cell_m(1, 3) = 0.0_dp
     368          10 :                helium_env(k)%helium%cell_m(2, 3) = 0.0_dp
     369          10 :                helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size
     370          10 :                helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size
     371          10 :                helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp
     372          10 :                helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp
     373          10 :                helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp
     374          10 :                helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size
     375          10 :                helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp
     376          10 :                helium_env(k)%helium%cell_m_inv(1, 3) = 0.0_dp
     377          10 :                helium_env(k)%helium%cell_m_inv(2, 3) = 0.0_dp
     378          10 :                helium_env(k)%helium%cell_m_inv(3, 3) = 1.0_dp/helium_env(k)%helium%cell_size
     379             :             CASE DEFAULT
     380           0 :                helium_env(k)%helium%cell_m(:, :) = 0.0_dp
     381          30 :                helium_env(k)%helium%cell_m_inv(:, :) = 0.0_dp
     382             : 
     383             :             END SELECT
     384             : 
     385             :          END DO ! k
     386             : 
     387          24 :          IF (logger%para_env%is_source()) THEN
     388             :             CALL section_vals_val_get(helium_section, "POTENTIAL_FILE_NAME", &
     389          12 :                                       c_val=potential_file_name)
     390             :             CALL open_file(file_name=TRIM(potential_file_name), &
     391          12 :                            file_action="READ", file_status="OLD", unit_number=input_unit)
     392          12 :             READ (input_unit, "(A)") msg_str
     393          12 :             READ (msg_str, *, IOSTAT=i) nlines, pdx, tau, &
     394          24 :                x_spline, dx, he_mass
     395          12 :             IF (i /= 0) THEN
     396             :                ! old style potential file, use default mass and potential
     397          12 :                he_mass = 4.00263037059764_dp !< 4He mass in [u]
     398          12 :                expl_pot = .FALSE.
     399          12 :                READ (msg_str, *, IOSTAT=i) nlines, pdx, tau, &
     400          24 :                   x_spline, dx
     401          12 :                IF (i /= 0) THEN
     402           0 :                   msg_str = "Format/Read Error from Solvent POTENTIAL_FILE"
     403           0 :                   CPABORT(msg_str)
     404             :                END IF
     405             :             ELSE
     406           0 :                expl_pot = .TRUE.
     407             :                ! in file really hb2m in kelvin angstrom**2
     408           0 :                he_mass = angstrom**2*kelvin/massunit/he_mass
     409             :                ! tau might be negative to get older versions of cp2k,
     410             :                ! that cannot handle the new potential file format to
     411             :                ! crash and not run a calculation with wrong mass/potential
     412           0 :                tau = ABS(tau)
     413             :             END IF
     414          12 :             tau = kelvin/tau
     415          12 :             x_spline = x_spline/angstrom
     416          12 :             dx = dx/angstrom
     417             :          END IF
     418             : 
     419          24 :          CALL helium_env(1)%comm%bcast(nlines, logger%para_env%source)
     420          24 :          CALL helium_env(1)%comm%bcast(pdx, logger%para_env%source)
     421          24 :          CALL helium_env(1)%comm%bcast(tau, logger%para_env%source)
     422          24 :          CALL helium_env(1)%comm%bcast(x_spline, logger%para_env%source)
     423          24 :          CALL helium_env(1)%comm%bcast(dx, logger%para_env%source)
     424          24 :          CALL helium_env(1)%comm%bcast(he_mass, logger%para_env%source)
     425          24 :          isize = (pdx + 1)*(pdx + 2) + 1
     426          96 :          ALLOCATE (pot_transfer(nlines, isize))
     427          24 :          IF (logger%para_env%is_source()) THEN
     428          12 :             IF (expl_pot) THEN
     429           0 :                DO i = 1, nlines
     430           0 :                   READ (input_unit, *) pot_transfer(i, :)
     431             :                END DO
     432             :             ELSE
     433        1212 :                DO i = 1, nlines
     434        1200 :                   READ (input_unit, *) pot_transfer(i, 1:isize - 1)
     435             :                   ! potential implicit, calculate it here now:
     436        1212 :                   pot_transfer(i, isize) = helium_vij(x_spline + (i - 1)*dx)*kelvin
     437             :                END DO
     438             :             END IF
     439          12 :             CALL close_file(unit_number=input_unit)
     440             :          END IF
     441          24 :          CALL helium_env(1)%comm%bcast(pot_transfer, logger%para_env%source)
     442             : 
     443          24 :          CALL spline_data_create(helium_env(1)%helium%vij)
     444          24 :          CALL init_splinexy(helium_env(1)%helium%vij, nlines)
     445          24 :          helium_env(1)%helium%vij%x1 = x_spline
     446             : 
     447          24 :          CALL spline_data_create(helium_env(1)%helium%u0)
     448          24 :          CALL init_splinexy(helium_env(1)%helium%u0, nlines)
     449          24 :          helium_env(1)%helium%u0%x1 = x_spline
     450             : 
     451          24 :          CALL spline_data_create(helium_env(1)%helium%e0)
     452          24 :          CALL init_splinexy(helium_env(1)%helium%e0, nlines)
     453          24 :          helium_env(1)%helium%e0%x1 = x_spline
     454             : 
     455          24 :          isize = pdx + 1
     456          24 :          ntab = ((isize + 1)*isize)/2 - 1   ! -1 because endpoint approx treated separately
     457          96 :          ALLOCATE (helium_env(1)%helium%uoffdiag(ntab, 2, nlines))
     458          96 :          ALLOCATE (helium_env(1)%helium%eoffdiag(ntab, 2, nlines))
     459          96 :          DO j = 1, isize
     460         240 :             DO i = j, isize
     461         144 :                IF (i + j == 2) CYCLE ! endpoint approx later separately
     462         120 :                k = ((i - 1)*i)/2 + j
     463       12120 :                helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)
     464         120 :                CALL init_spline(helium_env(1)%helium%vij, dx=dx)
     465       12120 :                helium_env(1)%helium%uoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:)
     466       12120 :                helium_env(1)%helium%uoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:)
     467         120 :                k = k + ((isize + 1)*isize)/2
     468       12120 :                helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)/kelvin
     469         120 :                CALL init_spline(helium_env(1)%helium%vij, dx=dx)
     470       12120 :                helium_env(1)%helium%eoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:)
     471       12120 :                helium_env(1)%helium%eoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:)
     472         216 :                ntab = ntab - 1
     473             :             END DO
     474             :          END DO
     475             : 
     476          24 :          ntab = SIZE(pot_transfer, 2)
     477        2424 :          helium_env(1)%helium%vij%y(:) = pot_transfer(:, ntab)/kelvin
     478          24 :          CALL init_spline(helium_env(1)%helium%vij, dx=dx)
     479             : 
     480        2424 :          helium_env(1)%helium%u0%y(:) = pot_transfer(:, 1)
     481          24 :          CALL init_spline(helium_env(1)%helium%u0, dx=dx)
     482          24 :          k = ((isize + 1)*isize)/2 + 1
     483        2424 :          helium_env(1)%helium%e0%y(:) = pot_transfer(:, k)/kelvin
     484          24 :          CALL init_spline(helium_env(1)%helium%e0, dx=dx)
     485             : 
     486          30 :          DO k = 2, mepos
     487           6 :             helium_env(k)%helium%vij => helium_env(1)%helium%vij
     488           6 :             helium_env(k)%helium%u0 => helium_env(1)%helium%u0
     489           6 :             helium_env(k)%helium%e0 => helium_env(1)%helium%e0
     490           6 :             helium_env(k)%helium%uoffdiag => helium_env(1)%helium%uoffdiag
     491          30 :             helium_env(k)%helium%eoffdiag => helium_env(1)%helium%eoffdiag
     492             :          END DO
     493             : 
     494          54 :          DO k = 1, mepos
     495             : 
     496          30 :             helium_env(k)%helium%pdx = pdx
     497          30 :             helium_env(k)%helium%tau = tau
     498             : 
     499             :             ! boltzmann : Boltzmann constant [J/K]
     500             :             ! h_bar     : Planck constant [J*s]
     501             :             ! J = kg*m^2/s^2
     502             :             ! 4He mass in [kg]
     503          30 :             mHe = he_mass*a_mass
     504             :             ! physical temperature [K]
     505          30 :             T = kelvin/helium_env(k)%helium%tau/helium_env(k)%helium%beads
     506             :             ! prefactors for calculating superfluid fractions [Angstrom^-2]
     507          30 :             helium_env(k)%helium%wpref = (((1e-20/h_bar)*mHe)/h_bar)*boltzmann*T
     508          30 :             helium_env(k)%helium%apref = (((4e-20/h_bar)*mHe)/h_bar)*boltzmann*T
     509             : 
     510          30 :             helium_env(k)%helium%he_mass_au = he_mass*massunit
     511          30 :             helium_env(k)%helium%hb2m = 1.0_dp/helium_env(k)%helium%he_mass_au
     512          30 :             helium_env(k)%helium%pweight = 0.0_dp
     513             : 
     514             :             ! Default in case sampling_method is not helium_sampling_worm.
     515          30 :             helium_env(k)%helium%worm_max_open_cycles = 0
     516             : 
     517             :             ! Choose sampling method:
     518             :             CALL section_vals_val_get(helium_section, "SAMPLING_METHOD", &
     519          30 :                                       i_val=helium_env(k)%helium%sampling_method)
     520             : 
     521          52 :             SELECT CASE (helium_env(k)%helium%sampling_method)
     522             :             CASE (helium_sampling_ceperley)
     523             :                ! check value of maxcycle
     524             :                CALL section_vals_val_get(helium_section, "CEPERLEY%MAX_PERM_CYCLE", &
     525          22 :                                          i_val=helium_env(k)%helium%maxcycle)
     526          22 :                i = helium_env(k)%helium%maxcycle
     527          22 :                CPASSERT(i >= 0)
     528          22 :                i = helium_env(k)%helium%atoms - helium_env(k)%helium%maxcycle
     529          22 :                CPASSERT(i >= 0)
     530             : 
     531             :                ! set m-distribution parameters
     532             :                CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%DISTRIBUTION-TYPE", &
     533          22 :                                          i_val=i)
     534          22 :                CPASSERT(i >= 1)
     535          22 :                CPASSERT(i <= 6)
     536          22 :                helium_env(k)%helium%m_dist_type = i
     537             :                CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-VALUE", &
     538          22 :                                          i_val=i)
     539          22 :                CPASSERT(i >= 1)
     540          22 :                CPASSERT(i <= helium_env(k)%helium%maxcycle)
     541          22 :                helium_env(k)%helium%m_value = i
     542             :                CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-RATIO", &
     543          22 :                                          r_val=rtmp)
     544          22 :                CPASSERT(rtmp > 0.0_dp)
     545          22 :                CPASSERT(rtmp <= 1.0_dp)
     546          22 :                helium_env(k)%helium%m_ratio = rtmp
     547             : 
     548             :                CALL section_vals_val_get(helium_section, "CEPERLEY%BISECTION", &
     549          22 :                                          i_val=helium_env(k)%helium%bisection)
     550             :                ! precheck bisection value (not all invalids are filtered out here yet)
     551          22 :                i = helium_env(k)%helium%bisection
     552          22 :                CPASSERT(i > 1)
     553          22 :                i = helium_env(k)%helium%beads - helium_env(k)%helium%bisection
     554          22 :                CPASSERT(i > 0)
     555             :                !
     556          22 :                itmp = helium_env(k)%helium%bisection
     557          22 :                rtmp = 2.0_dp**(ANINT(LOG(REAL(itmp, dp))/LOG(2.0_dp)))
     558          22 :                tcheck = ABS(REAL(itmp, KIND=dp) - rtmp)
     559          22 :                IF (tcheck > 100.0_dp*EPSILON(0.0_dp)) THEN
     560           0 :                   msg_str = "BISECTION should be integer power of 2."
     561           0 :                   CPABORT(msg_str)
     562             :                END IF
     563          22 :                helium_env(k)%helium%bisctlog2 = NINT(LOG(REAL(itmp, dp))/LOG(2.0_dp))
     564             : 
     565             :             CASE (helium_sampling_worm)
     566           8 :                NULLIFY (input_worm)
     567             :                input_worm => section_vals_get_subs_vals(helium_env(k)%helium%input, &
     568           8 :                                                         "MOTION%PINT%HELIUM%WORM")
     569             :                CALL section_vals_val_get(helium_section, "WORM%CENTROID_DRMAX", &
     570           8 :                                          r_val=helium_env(k)%helium%worm_centroid_drmax)
     571             : 
     572             :                CALL section_vals_val_get(helium_section, "WORM%STAGING_L", &
     573           8 :                                          i_val=helium_env(k)%helium%worm_staging_l)
     574             : 
     575             :                CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_SCALE", &
     576           8 :                                          r_val=helium_env(k)%helium%worm_open_close_scale)
     577             : 
     578             :                CALL section_vals_val_get(helium_section, "WORM%ALLOW_OPEN", &
     579           8 :                                          l_val=helium_env(k)%helium%worm_allow_open)
     580           8 :                IF (helium_env(k)%helium%atoms == 1) THEN
     581           0 :                   IF (helium_env(k)%helium%worm_allow_open) THEN
     582             :                      msg_str = "Default enabled open state sampling "// &
     583           0 :                                "for only 1 He might be inefficient."
     584           0 :                      CPWARN(msg_str)
     585             :                   END IF
     586             :                END IF
     587             : 
     588             :                CALL section_vals_val_get(helium_section, "WORM%MAX_OPEN_CYCLES", &
     589           8 :                                          i_val=helium_env(k)%helium%worm_max_open_cycles)
     590             : 
     591           8 :                IF (helium_env(k)%helium%worm_staging_l + 1 >= helium_env(k)%helium%beads) THEN
     592           0 :                   msg_str = "STAGING_L for worm sampling is too large"
     593           0 :                   CPABORT(msg_str)
     594           8 :                ELSE IF (helium_env(k)%helium%worm_staging_l < 1) THEN
     595           0 :                   msg_str = "STAGING_L must be positive integer"
     596           0 :                   CPABORT(msg_str)
     597             :                END IF
     598             : 
     599             :                CALL section_vals_val_get(helium_section, "WORM%SHOW_STATISTICS", &
     600           8 :                                          l_val=helium_env(k)%helium%worm_show_statistics)
     601             : 
     602             :                ! precompute an expensive scaling for the open and close acceptance probability
     603             :                ! tau is not included here, as It will be first defined in a few lines
     604           8 :                rtmp = 2.0_dp*pi*helium_env(k)%helium%hb2m
     605           8 :                rtmp = SQRT(rtmp)
     606           8 :                rtmp = rtmp**3
     607           8 :                rtmp = rtmp*helium_env(k)%helium%worm_open_close_scale
     608           8 :                IF (helium_env(k)%helium%periodic) THEN
     609           4 :                   rtmp = rtmp*helium_env(k)%helium%density
     610             :                ELSE
     611             :                   rtmp = rtmp*helium_env(k)%helium%atoms/ &
     612           4 :                          (4.0_dp/3.0_dp*pi*helium_env(k)%helium%droplet_radius**3)
     613             :                END IF
     614           8 :                helium_env(k)%helium%worm_ln_openclose_scale = LOG(rtmp)
     615             : 
     616             :                ! deal with acceptance statistics without changing the ceperley stuff
     617           8 :                helium_env(k)%helium%maxcycle = 1
     618           8 :                helium_env(k)%helium%bisctlog2 = 0
     619             : 
     620             :                ! get the absolute weights of the individual moves
     621           8 :                helium_env(k)%helium%worm_all_limit = 0
     622             :                CALL section_vals_val_get(helium_section, "WORM%CENTROID_WEIGHT", &
     623           8 :                                          i_val=itmp)
     624           8 :                helium_env(k)%helium%worm_centroid_min = 1
     625           8 :                helium_env(k)%helium%worm_centroid_max = itmp
     626           8 :                helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     627             :                CALL section_vals_val_get(helium_section, "WORM%STAGING_WEIGHT", &
     628           8 :                                          i_val=itmp)
     629           8 :                helium_env(k)%helium%worm_staging_min = helium_env(k)%helium%worm_centroid_max + 1
     630           8 :                helium_env(k)%helium%worm_staging_max = helium_env(k)%helium%worm_centroid_max + itmp
     631           8 :                helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     632           8 :                IF (helium_env(k)%helium%worm_allow_open) THEN
     633             :                   CALL section_vals_val_get(helium_section, "WORM%CRAWL_WEIGHT", &
     634           8 :                                             i_val=itmp)
     635           8 :                   helium_env(k)%helium%worm_fcrawl_min = helium_env(k)%helium%worm_staging_max + 1
     636           8 :                   helium_env(k)%helium%worm_fcrawl_max = helium_env(k)%helium%worm_staging_max + itmp
     637           8 :                   helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     638           8 :                   helium_env(k)%helium%worm_bcrawl_min = helium_env(k)%helium%worm_fcrawl_max + 1
     639           8 :                   helium_env(k)%helium%worm_bcrawl_max = helium_env(k)%helium%worm_fcrawl_max + itmp
     640           8 :                   helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     641             :                   CALL section_vals_val_get(helium_section, "WORM%HEAD_TAIL_WEIGHT", &
     642           8 :                                             i_val=itmp)
     643           8 :                   helium_env(k)%helium%worm_head_min = helium_env(k)%helium%worm_bcrawl_max + 1
     644           8 :                   helium_env(k)%helium%worm_head_max = helium_env(k)%helium%worm_bcrawl_max + itmp
     645           8 :                   helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     646           8 :                   helium_env(k)%helium%worm_tail_min = helium_env(k)%helium%worm_head_max + 1
     647           8 :                   helium_env(k)%helium%worm_tail_max = helium_env(k)%helium%worm_head_max + itmp
     648           8 :                   helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     649             :                   CALL section_vals_val_get(helium_section, "WORM%SWAP_WEIGHT", &
     650           8 :                                             i_val=itmp)
     651           8 :                   helium_env(k)%helium%worm_swap_min = helium_env(k)%helium%worm_tail_max + 1
     652           8 :                   helium_env(k)%helium%worm_swap_max = helium_env(k)%helium%worm_tail_max + itmp
     653           8 :                   helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     654             :                   CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_WEIGHT", &
     655           8 :                                             i_val=itmp)
     656           8 :                   helium_env(k)%helium%worm_open_close_min = helium_env(k)%helium%worm_swap_max + 1
     657           8 :                   helium_env(k)%helium%worm_open_close_max = helium_env(k)%helium%worm_swap_max + itmp
     658           8 :                   helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
     659             :                   CALL section_vals_val_get(helium_section, "WORM%CRAWL_REPETITION", &
     660           8 :                                             i_val=helium_env(k)%helium%worm_repeat_crawl)
     661             :                END IF
     662             : 
     663             :                !CPPostcondition(i<helium_env(k)%helium%beads,cp_failure_level,routineP,failure)
     664             :                ! end of worm
     665             :             CASE DEFAULT
     666           0 :                msg_str = "Unknown helium sampling method!"
     667         104 :                CPABORT(msg_str)
     668             :             END SELECT
     669             : 
     670             :             ! ALLOCATE helium-related arrays
     671          30 :             i = helium_env(k)%helium%atoms
     672          30 :             j = helium_env(k)%helium%beads
     673         120 :             ALLOCATE (helium_env(k)%helium%pos(3, i, j))
     674       56280 :             helium_env(k)%helium%pos = 0.0_dp
     675         120 :             ALLOCATE (helium_env(k)%helium%work(3, i, j))
     676          90 :             ALLOCATE (helium_env(k)%helium%ptable(helium_env(k)%helium%maxcycle + 1))
     677          90 :             ALLOCATE (helium_env(k)%helium%permutation(i))
     678          90 :             ALLOCATE (helium_env(k)%helium%iperm(i))
     679         120 :             ALLOCATE (helium_env(k)%helium%tmatrix(i, i))
     680         120 :             ALLOCATE (helium_env(k)%helium%nmatrix(i, 2*i))
     681         120 :             ALLOCATE (helium_env(k)%helium%pmatrix(i, i))
     682         120 :             ALLOCATE (helium_env(k)%helium%ipmatrix(i, i))
     683          30 :             itmp = helium_env(k)%helium%bisctlog2 + 2
     684         120 :             ALLOCATE (helium_env(k)%helium%num_accepted(itmp, helium_env(k)%helium%maxcycle))
     685          90 :             ALLOCATE (helium_env(k)%helium%plength_avrg(helium_env(k)%helium%atoms))
     686          90 :             ALLOCATE (helium_env(k)%helium%plength_inst(helium_env(k)%helium%atoms))
     687          90 :             ALLOCATE (helium_env(k)%helium%atom_plength(helium_env(k)%helium%atoms))
     688          30 :             IF (helium_env(k)%helium%worm_max_open_cycles > 0) THEN
     689           0 :                ALLOCATE (helium_env(k)%helium%savepermutation(i))
     690           0 :                ALLOCATE (helium_env(k)%helium%saveiperm(i))
     691           0 :                ALLOCATE (helium_env(k)%helium%savepos(3, i, j))
     692             :             END IF
     693             : 
     694             :             ! check whether rdfs should be calculated and printed
     695          30 :             helium_env(k)%helium%rdf_present = helium_property_active(helium_env(k)%helium, "RDF")
     696          30 :             IF (helium_env(k)%helium%rdf_present) THEN
     697             :                ! allocate & initialize rdf related data structures
     698          30 :                CALL helium_rdf_init(helium_env(k)%helium)
     699             :             END IF
     700             : 
     701             :             ! check whether densities should be calculated and printed
     702          30 :             helium_env(k)%helium%rho_present = helium_property_active(helium_env(k)%helium, "RHO")
     703          54 :             IF (helium_env(k)%helium%rho_present) THEN
     704             :                ! allocate & initialize density related data structures
     705          30 :                NULLIFY (helium_env(k)%helium%rho_property)
     706          30 :                CALL helium_rho_init(helium_env(k)%helium)
     707             :             END IF
     708             : 
     709             :          END DO
     710             : 
     711             :          ! restore averages calculated in previous runs
     712          24 :          CALL helium_averages_restore(helium_env)
     713             : 
     714          54 :          DO k = 1, mepos
     715             :             ! fill in the solute-related data structures
     716          30 :             helium_env(k)%helium%e_corr = 0.0_dp
     717          30 :             IF (helium_env(k)%helium%solute_present) THEN
     718          20 :                IF (helium_env(k)%helium%solute_beads > helium_env(k)%helium%beads) THEN
     719             :                   ! Imaginary time striding for solute:
     720             :                   helium_env(k)%helium%bead_ratio = helium_env(k)%helium%solute_beads/ &
     721           2 :                                                     helium_env(k)%helium%beads
     722             :                   ! check if bead numbers are commensurate:
     723           2 :                   i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%beads - helium_env(k)%helium%solute_beads
     724           2 :                   IF (i /= 0) THEN
     725           0 :                      msg_str = "Adjust number of solute beads to multiple of solvent beads."
     726           0 :                      CPABORT(msg_str)
     727             :                   END IF
     728             :                   msg_str = "Using multiple-time stepping in imaginary time for solute to couple "// &
     729             :                             "to fewer solvent beads, e.g. for factor 3: "// &
     730             :                             "he_1 - 3*sol_1; he_2 - 3*sol_4... "// &
     731           2 :                             "Avoid too large coupling factors."
     732           2 :                   CPWARN(msg_str)
     733          18 :                ELSE IF (helium_env(k)%helium%solute_beads < helium_env(k)%helium%beads) THEN
     734             :                   ! Imaginary time striding for solvent:
     735             :                   helium_env(k)%helium%bead_ratio = helium_env(k)%helium%beads/ &
     736          18 :                                                     helium_env(k)%helium%solute_beads
     737             :                   ! check if bead numbers are commensurate:
     738          18 :                   i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%solute_beads - helium_env(k)%helium%beads
     739          18 :                   IF (i /= 0) THEN
     740           0 :                      msg_str = "Adjust number of solvent beads to multiple of solute beads."
     741           0 :                      CPABORT(msg_str)
     742             :                   END IF
     743             :                   msg_str = "Coupling solvent beads to fewer solute beads via "// &
     744             :                             "direct coupling, e.g. for factor 3: "// &
     745          18 :                             "sol_1 - he_1,2,3; sol_2 - he_4,5,6..."
     746          18 :                   CPWARN(msg_str)
     747             :                END IF
     748             : !TODO       Adjust helium bead number if not comm. and if coords not given expl.
     749             : 
     750             :                ! check if tau, temperature and bead number are consistent:
     751          20 :                tcheck = ABS((helium_env(k)%helium%tau*helium_env(k)%helium%beads - solute%beta)/solute%beta)
     752          20 :                IF (tcheck > 1.0e-14_dp) THEN
     753           0 :                   msg_str = "Tau, temperature and bead number are inconsistent."
     754           0 :                   CPABORT(msg_str)
     755             :                END IF
     756             : 
     757          20 :                CALL helium_set_solute_indices(helium_env(k)%helium, solute)
     758          20 :                CALL helium_set_solute_cell(helium_env(k)%helium, solute)
     759             : 
     760             :                ! set the interaction potential type
     761             :                CALL section_vals_val_get(helium_section, "SOLUTE_INTERACTION", &
     762          20 :                                          i_val=helium_env(k)%helium%solute_interaction)
     763          20 :                IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_none) THEN
     764             :                   WRITE (msg_str, '(A,I0,A)') &
     765             :                      "Solute found but no helium-solute interaction selected "// &
     766           0 :                      "(see SOLUTE_INTERACTION keyword)"
     767           0 :                   CPABORT(msg_str)
     768             :                END IF
     769             : 
     770             :                ! ALLOCATE solute-related arrays
     771             :                ALLOCATE (helium_env(k)%helium%force_avrg(helium_env(k)%helium%solute_beads, &
     772          80 :                                                          helium_env(k)%helium%solute_atoms*3))
     773             :                ALLOCATE (helium_env(k)%helium%force_inst(helium_env(k)%helium%solute_beads, &
     774          80 :                                                          helium_env(k)%helium%solute_atoms*3))
     775             : 
     776          60 :                ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_1d(solute%p*solute%ndim))
     777          60 :                ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_np_1d(solute%p*solute%ndim*helium_env(k)%helium%num_env))
     778          80 :                ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_2d(solute%p, solute%ndim))
     779             : 
     780             :             ELSE
     781          10 :                helium_env(k)%helium%bead_ratio = 0
     782          10 :                IF (helium_env(k)%helium%periodic) THEN
     783             :                   ! this assumes a specific potential (and its ugly):
     784          10 :                   x1 = angstrom*0.5_dp*helium_env(k)%helium%cell_size
     785             :                   ! 10.8 is in Kelvin, x1 needs to be in Angstrom,
     786             :                   ! since 2.9673 is in Angstrom
     787             :                   helium_env(k)%helium%e_corr = (twopi*helium_env(k)%helium%density/angstrom**3*10.8_dp* &
     788             :                                                  (544850.4_dp*EXP(-13.353384_dp*x1/2.9673_dp)*(2.9673_dp/13.353384_dp)**3* &
     789             :                                                   (2.0_dp + 2.0_dp*13.353384_dp*x1/2.9673_dp + (13.353384_dp*x1/2.9673_dp)**2) - &
     790             :                                                   (((0.1781_dp/7.0_dp*(2.9673_dp/x1)**2 + 0.4253785_dp/5.0_dp)*(2.9673_dp/x1)**2 + &
     791          10 :                                                     1.3732412_dp/3.0_dp)*(2.9673_dp/x1)**3)*2.9673_dp**3))/kelvin
     792             :                END IF
     793             :             END IF
     794             : 
     795             :             ! ALLOCATE temporary arrays
     796          90 :             ALLOCATE (helium_env(k)%helium%itmp_atoms_1d(helium_env(k)%helium%atoms))
     797          90 :             ALLOCATE (helium_env(k)%helium%ltmp_atoms_1d(helium_env(k)%helium%atoms))
     798          90 :             ALLOCATE (helium_env(k)%helium%itmp_atoms_np_1d(helium_env(k)%helium%atoms*helium_env(k)%helium%num_env))
     799          90 :             ALLOCATE (helium_env(k)%helium%rtmp_3_np_1d(3*helium_env(k)%helium%num_env))
     800             :             ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_1d(3*helium_env(k)%helium%atoms* &
     801          90 :                                                                  helium_env(k)%helium%beads))
     802             :             ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_np_1d(3*helium_env(k)%helium%atoms* &
     803          90 :                                                                     helium_env(k)%helium%beads*helium_env(k)%helium%num_env))
     804             :             ALLOCATE (helium_env(k)%helium%ltmp_3_atoms_beads_3d(3, helium_env(k)%helium%atoms, &
     805         120 :                                                                  helium_env(k)%helium%beads))
     806          54 :             IF (k .EQ. 1) THEN
     807          24 :                CALL helium_write_setup(helium_env(k)%helium)
     808             :             END IF
     809             :          END DO
     810          24 :          DEALLOCATE (pot_transfer)
     811             :       ELSE
     812             :          ! Deallocate env_all on processors without helium_env
     813           0 :          DEALLOCATE (env_all)
     814             :       END IF ! mepos .GT. 0
     815             : 
     816          24 :       NULLIFY (env_all)
     817          24 :       CALL timestop(handle)
     818             : 
     819          24 :       RETURN
     820          96 :    END SUBROUTINE helium_create
     821             : 
     822             : ! ***************************************************************************
     823             : !> \brief  Releases helium_solvent_type
     824             : !> \param helium_env ...
     825             : !> \par    History
     826             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     827             : !> \author hforbert
     828             : ! **************************************************************************************************
     829          24 :    SUBROUTINE helium_release(helium_env)
     830             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     831             : 
     832             :       INTEGER                                            :: k
     833             : 
     834          24 :       IF (ASSOCIATED(helium_env)) THEN
     835          54 :          DO k = 1, SIZE(helium_env)
     836          30 :             IF (k .EQ. 1) THEN
     837          24 :                CALL helium_env(k)%comm%free()
     838          24 :                DEALLOCATE (helium_env(k)%env_all)
     839             :             END IF
     840          30 :             NULLIFY (helium_env(k)%env_all)
     841             : 
     842             :             ! DEALLOCATE temporary arrays
     843           0 :             DEALLOCATE ( &
     844             :                helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
     845           0 :                helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
     846           0 :                helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
     847           0 :                helium_env(k)%helium%rtmp_3_np_1d, &
     848           0 :                helium_env(k)%helium%itmp_atoms_np_1d, &
     849           0 :                helium_env(k)%helium%ltmp_atoms_1d, &
     850          30 :                helium_env(k)%helium%itmp_atoms_1d)
     851             : 
     852             :             NULLIFY ( &
     853             :                helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
     854          30 :                helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
     855          30 :                helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
     856          30 :                helium_env(k)%helium%rtmp_3_np_1d, &
     857          30 :                helium_env(k)%helium%itmp_atoms_np_1d, &
     858          30 :                helium_env(k)%helium%ltmp_atoms_1d, &
     859          30 :                helium_env(k)%helium%itmp_atoms_1d &
     860          30 :                )
     861             : 
     862          30 :             IF (helium_env(k)%helium%solute_present) THEN
     863             :                ! DEALLOCATE solute-related arrays
     864           0 :                DEALLOCATE ( &
     865             :                   helium_env(k)%helium%rtmp_p_ndim_2d, &
     866           0 :                   helium_env(k)%helium%rtmp_p_ndim_np_1d, &
     867           0 :                   helium_env(k)%helium%rtmp_p_ndim_1d, &
     868           0 :                   helium_env(k)%helium%force_inst, &
     869          20 :                   helium_env(k)%helium%force_avrg)
     870             :                NULLIFY ( &
     871             :                   helium_env(k)%helium%rtmp_p_ndim_2d, &
     872          20 :                   helium_env(k)%helium%rtmp_p_ndim_np_1d, &
     873          20 :                   helium_env(k)%helium%rtmp_p_ndim_1d, &
     874          20 :                   helium_env(k)%helium%force_inst, &
     875          20 :                   helium_env(k)%helium%force_avrg)
     876             :             END IF
     877             : 
     878          30 :             IF (helium_env(k)%helium%rho_present) THEN
     879           0 :                DEALLOCATE ( &
     880             :                   helium_env(k)%helium%rho_rstr, &
     881           0 :                   helium_env(k)%helium%rho_accu, &
     882           0 :                   helium_env(k)%helium%rho_inst, &
     883          30 :                   helium_env(k)%helium%rho_incr)
     884             :                NULLIFY ( &
     885             :                   helium_env(k)%helium%rho_rstr, &
     886          30 :                   helium_env(k)%helium%rho_accu, &
     887          30 :                   helium_env(k)%helium%rho_inst, &
     888          30 :                   helium_env(k)%helium%rho_incr)
     889             :                ! DEALLOCATE everything
     890          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix)
     891          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_name)
     892          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_index)
     893          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix)
     894          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_name)
     895          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_index)
     896          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix)
     897          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_name)
     898          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_index)
     899          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix)
     900          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_name)
     901          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_index)
     902          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix)
     903          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name)
     904          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index)
     905          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix)
     906          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name)
     907          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index)
     908          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix)
     909          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_name)
     910          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_index)
     911          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix)
     912          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_name)
     913          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_index)
     914          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix)
     915          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name)
     916          30 :                DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index)
     917          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix)
     918          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name)
     919          30 :                NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index)
     920          30 :                DEALLOCATE (helium_env(k)%helium%rho_property)
     921          30 :                NULLIFY (helium_env(k)%helium%rho_property)
     922             :             END IF
     923             : 
     924          30 :             CALL helium_rdf_release(helium_env(k)%helium)
     925             : 
     926             :             ! DEALLOCATE helium-related arrays
     927           0 :             DEALLOCATE ( &
     928             :                helium_env(k)%helium%atom_plength, &
     929           0 :                helium_env(k)%helium%plength_inst, &
     930           0 :                helium_env(k)%helium%plength_avrg, &
     931           0 :                helium_env(k)%helium%num_accepted, &
     932           0 :                helium_env(k)%helium%ipmatrix, &
     933           0 :                helium_env(k)%helium%pmatrix, &
     934           0 :                helium_env(k)%helium%nmatrix, &
     935           0 :                helium_env(k)%helium%tmatrix, &
     936           0 :                helium_env(k)%helium%iperm, &
     937           0 :                helium_env(k)%helium%permutation, &
     938           0 :                helium_env(k)%helium%ptable, &
     939           0 :                helium_env(k)%helium%work, &
     940          30 :                helium_env(k)%helium%pos)
     941          30 :             IF (helium_env(k)%helium%worm_max_open_cycles > 0) THEN
     942           0 :                DEALLOCATE (helium_env(k)%helium%saveiperm, &
     943           0 :                            helium_env(k)%helium%savepermutation, &
     944           0 :                            helium_env(k)%helium%savepos)
     945             :             END IF
     946             :             NULLIFY ( &
     947             :                helium_env(k)%helium%atom_plength, &
     948          30 :                helium_env(k)%helium%plength_inst, &
     949          30 :                helium_env(k)%helium%plength_avrg, &
     950          30 :                helium_env(k)%helium%num_accepted, &
     951          30 :                helium_env(k)%helium%ipmatrix, &
     952          30 :                helium_env(k)%helium%pmatrix, &
     953          30 :                helium_env(k)%helium%nmatrix, &
     954          30 :                helium_env(k)%helium%tmatrix, &
     955          30 :                helium_env(k)%helium%iperm, &
     956          30 :                helium_env(k)%helium%saveiperm, &
     957          30 :                helium_env(k)%helium%permutation, &
     958          30 :                helium_env(k)%helium%savepermutation, &
     959          30 :                helium_env(k)%helium%ptable, &
     960          30 :                helium_env(k)%helium%work, &
     961          30 :                helium_env(k)%helium%pos, &
     962          30 :                helium_env(k)%helium%savepos &
     963          30 :                )
     964             : 
     965          30 :             IF (k == 1) THEN
     966          24 :                CALL spline_data_release(helium_env(k)%helium%vij)
     967          24 :                CALL spline_data_release(helium_env(k)%helium%u0)
     968          24 :                CALL spline_data_release(helium_env(k)%helium%e0)
     969          24 :                DEALLOCATE (helium_env(k)%helium%uoffdiag)
     970          24 :                DEALLOCATE (helium_env(k)%helium%eoffdiag)
     971             :             END IF
     972             :             NULLIFY (helium_env(k)%helium%uoffdiag, &
     973          30 :                      helium_env(k)%helium%eoffdiag, &
     974          30 :                      helium_env(k)%helium%vij, &
     975          30 :                      helium_env(k)%helium%u0, &
     976          30 :                      helium_env(k)%helium%e0)
     977             : 
     978          30 :             DEALLOCATE (helium_env(k)%helium%rng_stream_uniform)
     979          30 :             DEALLOCATE (helium_env(k)%helium%rng_stream_gaussian)
     980             : 
     981             :             ! deallocate solute-related arrays
     982          30 :             IF (helium_env(k)%helium%solute_present) THEN
     983          20 :                DEALLOCATE (helium_env(k)%helium%solute_element)
     984          20 :                NULLIFY (helium_env(k)%helium%solute_element)
     985             :             END IF
     986             : 
     987             :             ! Deallocate everything from the helium_set_solute_indices
     988          30 :             IF (ASSOCIATED(helium_env(k)%helium%ename)) THEN
     989           0 :                DEALLOCATE (helium_env(k)%helium%ename)
     990           0 :                NULLIFY (helium_env(k)%helium%ename)
     991             :             END IF
     992             : 
     993          54 :             DEALLOCATE (helium_env(k)%helium)
     994             : 
     995             :          END DO
     996             : 
     997          24 :          DEALLOCATE (helium_env)
     998             :          NULLIFY (helium_env)
     999             :       END IF
    1000          24 :       RETURN
    1001             :    END SUBROUTINE helium_release
    1002             : 
    1003             : ! ***************************************************************************
    1004             : !> \brief  Initialize helium data structures.
    1005             : !> \param helium_env ...
    1006             : !> \param pint_env ...
    1007             : !> \par    History
    1008             : !>         removed references to pint_env_type data structure [lwalewski]
    1009             : !>         2009-11-10 init/restore coords, perm, RNG and forces [lwalewski]
    1010             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1011             : !> \author hforbert
    1012             : !> \note   Initializes helium coordinates either as random positions or from
    1013             : !>         HELIUM%COORD section if it's present in the input file.
    1014             : !>         Initializes helium permutation state as identity permutation or
    1015             : !>         from HELIUM%PERM section if it's present in the input file.
    1016             : ! **************************************************************************************************
    1017          24 :    SUBROUTINE helium_init(helium_env, pint_env)
    1018             : 
    1019             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1020             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    1021             : 
    1022             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_init'
    1023             : 
    1024             :       INTEGER                                            :: handle, k
    1025             :       LOGICAL                                            :: coords_presampled, explicit, presample
    1026             :       REAL(KIND=dp)                                      :: initkT, solute_radius
    1027             :       TYPE(cp_logger_type), POINTER                      :: logger
    1028             :       TYPE(section_vals_type), POINTER                   :: helium_section, sec
    1029             : 
    1030          24 :       CALL timeset(routineN, handle)
    1031             : 
    1032          24 :       NULLIFY (logger)
    1033          24 :       logger => cp_get_default_logger()
    1034             : 
    1035          24 :       IF (ASSOCIATED(helium_env)) THEN
    1036             : 
    1037          24 :          NULLIFY (helium_section)
    1038             :          helium_section => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1039          24 :                                                       "MOTION%PINT%HELIUM")
    1040             : 
    1041             :          ! restore RNG state
    1042          24 :          NULLIFY (sec)
    1043          24 :          sec => section_vals_get_subs_vals(helium_section, "RNG_STATE")
    1044          24 :          CALL section_vals_get(sec, explicit=explicit)
    1045          24 :          IF (explicit) THEN
    1046           6 :             CALL helium_rng_restore(helium_env)
    1047             :          ELSE
    1048          18 :             CALL helium_write_line("RNG state initialized as new.")
    1049             :          END IF
    1050             : 
    1051             :          ! init/restore permutation state
    1052          24 :          NULLIFY (sec)
    1053          24 :          sec => section_vals_get_subs_vals(helium_section, "PERM")
    1054          24 :          CALL section_vals_get(sec, explicit=explicit)
    1055          24 :          IF (explicit) THEN
    1056           6 :             CALL helium_perm_restore(helium_env)
    1057             :          ELSE
    1058          18 :             CALL helium_perm_init(helium_env)
    1059          18 :             CALL helium_write_line("Permutation state initialized as identity.")
    1060             :          END IF
    1061             : 
    1062             :          ! Specify if forces should be obtained as AVG or LAST
    1063          54 :          DO k = 1, SIZE(helium_env)
    1064             :             CALL section_vals_val_get(helium_section, "GET_FORCES", &
    1065          54 :                                       i_val=helium_env(k)%helium%get_helium_forces)
    1066             :          END DO
    1067             : 
    1068          54 :          DO k = 1, SIZE(helium_env)
    1069             :             ! init center of mass
    1070          54 :             IF (helium_env(k)%helium%solute_present) THEN
    1071          80 :                helium_env(k)%helium%center(:) = pint_com_pos(pint_env)
    1072             :             ELSE
    1073          40 :                helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
    1074             :             END IF
    1075             :          END DO
    1076             : 
    1077             :          ! init/restore coordinates
    1078          24 :          NULLIFY (sec)
    1079          24 :          sec => section_vals_get_subs_vals(helium_section, "COORD")
    1080          24 :          CALL section_vals_get(sec, explicit=explicit)
    1081          24 :          IF (explicit) THEN
    1082           6 :             CALL helium_coord_restore(helium_env)
    1083           6 :             CALL helium_write_line("Coordinates restarted.")
    1084             :          ELSE
    1085          18 :             CALL section_vals_val_get(helium_section, "COORD_INIT_TEMP", r_val=initkT)
    1086          18 :             CALL section_vals_val_get(helium_section, "SOLUTE_RADIUS", r_val=solute_radius)
    1087          18 :             CALL helium_coord_init(helium_env, initkT, solute_radius)
    1088          18 :             IF (initkT > 0.0_dp) THEN
    1089          12 :                CALL helium_write_line("Coordinates initialized with thermal gaussian.")
    1090             :             ELSE
    1091           6 :                CALL helium_write_line("Coordinates initialized as point particles.")
    1092             :             END IF
    1093             :          END IF
    1094             : 
    1095          54 :          DO k = 1, SIZE(helium_env)
    1096             : 
    1097          30 :             helium_env(k)%helium%worm_is_closed = .TRUE.
    1098          30 :             helium_env(k)%helium%worm_atom_idx = 0
    1099          30 :             helium_env(k)%helium%worm_bead_idx = 0
    1100             : 
    1101       56280 :             helium_env(k)%helium%work(:, :, :) = helium_env(k)%helium%pos(:, :, :)
    1102             : 
    1103             :             ! init center of mass
    1104          54 :             IF (helium_env(k)%helium%solute_present) THEN
    1105          80 :                helium_env(k)%helium%center(:) = pint_com_pos(pint_env)
    1106             :             ELSE
    1107          10 :                IF (helium_env(k)%helium%periodic) THEN
    1108          40 :                   helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
    1109             :                ELSE
    1110           0 :                   helium_env(k)%helium%center(:) = helium_com(helium_env(k)%helium)
    1111             :                END IF
    1112             :             END IF
    1113             :          END DO
    1114             : 
    1115             :          ! Optional helium coordinate presampling:
    1116             :          ! Assume IONODE to have always at least one helium_env
    1117             :          CALL section_vals_val_get(helium_section, "PRESAMPLE", &
    1118          24 :                                    l_val=presample)
    1119          24 :          coords_presampled = .FALSE.
    1120          24 :          IF (presample) THEN
    1121          30 :             DO k = 1, SIZE(helium_env)
    1122          30 :                helium_env(k)%helium%current_step = 0
    1123             :             END DO
    1124          12 :             CALL helium_sample(helium_env, pint_env)
    1125          30 :             DO k = 1, SIZE(helium_env)
    1126         558 :                IF (helium_env(k)%helium%solute_present) helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
    1127         198 :                helium_env(k)%helium%energy_avrg(:) = 0.0_dp
    1128         324 :                helium_env(k)%helium%plength_avrg(:) = 0.0_dp
    1129         228 :                helium_env(k)%helium%num_accepted(:, :) = 0.0_dp
    1130             :                ! Reset properties accumulated over presample:
    1131          72 :                helium_env(k)%helium%proarea%accu(:) = 0.0_dp
    1132          72 :                helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
    1133          72 :                helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
    1134          72 :                helium_env(k)%helium%mominer%accu(:) = 0.0_dp
    1135          18 :                IF (helium_env(k)%helium%rho_present) THEN
    1136    32165838 :                   helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
    1137             :                END IF
    1138          30 :                IF (helium_env(k)%helium%rdf_present) THEN
    1139       13014 :                   helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
    1140             :                END IF
    1141             :             END DO
    1142          12 :             coords_presampled = .TRUE.
    1143          12 :             CALL helium_write_line("Bead coordinates pre-sampled.")
    1144             :          END IF
    1145             : 
    1146          24 :          IF (helium_env(1)%helium%solute_present) THEN
    1147             :             ! restore helium forces
    1148          14 :             NULLIFY (sec)
    1149          14 :             sec => section_vals_get_subs_vals(helium_section, "FORCE")
    1150          14 :             CALL section_vals_get(sec, explicit=explicit)
    1151          14 :             IF (explicit) THEN
    1152           2 :                IF (.NOT. coords_presampled) THEN
    1153           2 :                   CALL helium_force_restore(helium_env)
    1154             :                END IF
    1155             :             ELSE
    1156          12 :                IF (.NOT. coords_presampled) THEN
    1157           6 :                   CALL helium_force_init(helium_env)
    1158           6 :                   CALL helium_write_line("Forces on the solute initialized as zero.")
    1159             :                END IF
    1160             :             END IF
    1161             :             !! Update pint_env force, assume always one helium_env at IONODE
    1162             :             !IF (pint_env%logger%para_env%is_source()) THEN
    1163             :             !   pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
    1164             :             !END IF
    1165             :             !CALL pint_env%logger%para_env%bcast(pint_env%f,&
    1166             :             !              pint_env%logger%para_env%source)
    1167             : 
    1168             :          END IF
    1169             :       END IF
    1170             : 
    1171          24 :       CALL timestop(handle)
    1172             : 
    1173          24 :       RETURN
    1174             :    END SUBROUTINE helium_init
    1175             : 
    1176             : ! ***************************************************************************
    1177             : ! Data transfer functions.
    1178             : !
    1179             : ! These functions manipulate and transfer data between the runtime
    1180             : ! environment and the input structure.
    1181             : ! ***************************************************************************
    1182             : 
    1183             : ! ***************************************************************************
    1184             : !> \brief  Initialize helium coordinates with random positions.
    1185             : !> \param helium_env ...
    1186             : !> \param initkT ...
    1187             : !> \param solute_radius ...
    1188             : !> \date   2009-11-09
    1189             : !> \par    History
    1190             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1191             : !>         2018-04-30 Useful initialization for droplets [fuhl]
    1192             : !> \author Lukasz Walewski
    1193             : ! **************************************************************************************************
    1194          18 :    SUBROUTINE helium_coord_init(helium_env, initkT, solute_radius)
    1195             :       TYPE(helium_solvent_p_type), DIMENSION(:), &
    1196             :          INTENT(INOUT), POINTER                          :: helium_env
    1197             :       REAL(KIND=dp), INTENT(IN)                          :: initkT, solute_radius
    1198             : 
    1199             :       REAL(KIND=dp), PARAMETER                           :: minHeHedst = 5.669177966_dp
    1200             : 
    1201             :       INTEGER                                            :: ia, ib, ic, id, iter, k
    1202             :       LOGICAL                                            :: invalidpos
    1203             :       REAL(KIND=dp)                                      :: minHeHedsttmp, r1, r2
    1204          18 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: centroids
    1205             :       REAL(KIND=dp), DIMENSION(3)                        :: cvek, rvek, tvek
    1206             : 
    1207             :       !corresponds to three angstrom (roughly first maximum of He-He-rdf)
    1208          18 :       minHeHedsttmp = minHeHedst
    1209             : 
    1210          42 :       DO k = 1, SIZE(helium_env)
    1211          24 :          IF (helium_env(k)%helium%solute_present) THEN
    1212          72 :             cvek(:) = helium_env(k)%helium%center(:)
    1213          18 :             CALL helium_pbc(helium_env(k)%helium, cvek)
    1214             :          END IF
    1215             : 
    1216          72 :          ALLOCATE (centroids(3, helium_env(k)%helium%atoms))
    1217          24 :          IF (helium_env(k)%helium%periodic) THEN
    1218         330 :             DO ia = 1, helium_env(k)%helium%atoms
    1219             :                invalidpos = .TRUE.
    1220             :                iter = 0
    1221      343843 :                DO WHILE (invalidpos)
    1222      343513 :                   iter = iter + 1
    1223      343513 :                   invalidpos = .FALSE.
    1224             :                   ! if sampling fails to often, reduce he he criterion
    1225             :                   !CS TODO:
    1226             :                   !minHeHedsttmp = 0.90_dp**(iter/100)*minHeHedst
    1227      343513 :                   minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst
    1228     1374052 :                   DO ic = 1, 3
    1229     1030539 :                      r1 = helium_env(k)%helium%rng_stream_uniform%next()
    1230     1030539 :                      r1 = 2.0_dp*r1 - 1.0_dp
    1231     1030539 :                      r1 = r1*helium_env(k)%helium%cell_size
    1232     1374052 :                      centroids(ic, ia) = r1
    1233             :                   END DO
    1234             :                   ! check if helium is outside of cell
    1235     1374052 :                   tvek(:) = centroids(:, ia)
    1236      343513 :                   CALL helium_pbc(helium_env(k)%helium, tvek(:))
    1237     1374052 :                   rvek(:) = tvek(:) - centroids(:, ia)
    1238     1374052 :                   r2 = DOT_PRODUCT(rvek, rvek)
    1239      343513 :                   IF (r2 > 1.0_dp*10.0_dp**(-6)) THEN
    1240             :                      invalidpos = .TRUE.
    1241             :                   ELSE
    1242             :                      ! check for helium-helium collision
    1243      201908 :                      DO id = 1, ia - 1
    1244      806348 :                         rvek = centroids(:, ia) - centroids(:, id)
    1245      201587 :                         CALL helium_pbc(helium_env(k)%helium, rvek)
    1246      806348 :                         r2 = DOT_PRODUCT(rvek, rvek)
    1247      201908 :                         IF (r2 < minHeHedsttmp**2) THEN
    1248             :                            invalidpos = .TRUE.
    1249             :                            EXIT
    1250             :                         END IF
    1251             :                      END DO
    1252             :                   END IF
    1253       21775 :                   IF (.NOT. invalidpos) THEN
    1254             :                      ! check if centroid collides with molecule
    1255         321 :                      IF (helium_env(k)%helium%solute_present) THEN
    1256         516 :                         rvek(:) = (cvek(:) - centroids(:, ia))
    1257         516 :                         r2 = DOT_PRODUCT(rvek, rvek)
    1258         129 :                         IF (r2 <= solute_radius**2) invalidpos = .TRUE.
    1259             :                      END IF
    1260             :                   END IF
    1261             :                END DO
    1262             :             END DO
    1263             :             ! do thermal gaussian delocalization of hot start
    1264          10 :             IF (initkT > 0.0_dp) THEN
    1265          10 :                CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkT)
    1266             :             ELSE
    1267           0 :                DO ia = 1, helium_env(k)%helium%atoms
    1268           0 :                   DO ib = 1, helium_env(k)%helium%beads
    1269           0 :                      helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia)
    1270             :                   END DO
    1271             :                END DO
    1272             :             END IF
    1273             :             ! apply PBC to bead coords
    1274         330 :             DO ia = 1, helium_env(k)%helium%atoms
    1275        7178 :                DO ib = 1, helium_env(k)%helium%beads
    1276        6848 :                   CALL helium_pbc(helium_env(k)%helium, helium_env(k)%helium%pos(:, ia, ib))
    1277             :                   ! check if bead collides with molecule
    1278        7168 :                   IF (helium_env(k)%helium%solute_present) THEN
    1279        8192 :                      rvek(:) = (cvek(:) - helium_env(k)%helium%pos(:, ia, ib))
    1280        8192 :                      r2 = DOT_PRODUCT(rvek, rvek)
    1281        2048 :                      IF (r2 <= solute_radius**2) THEN
    1282           2 :                         r1 = SQRT(r2)
    1283             :                         helium_env(k)%helium%pos(:, ia, ib) = &
    1284           8 :                            cvek(:) + solute_radius/r1*rvek(:)
    1285             :                      END IF
    1286             :                   END IF
    1287             :                END DO
    1288             :             END DO
    1289             :          ELSE
    1290         192 :             DO ia = 1, helium_env(k)%helium%atoms
    1291             :                iter = 0
    1292             :                invalidpos = .TRUE.
    1293         546 :                DO WHILE (invalidpos)
    1294         354 :                   invalidpos = .FALSE.
    1295         354 :                   iter = iter + 1
    1296             :                   ! if sampling fails to often, reduce he he criterion
    1297         354 :                   minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst
    1298        1416 :                   DO ic = 1, 3
    1299        1062 :                      rvek(ic) = helium_env(k)%helium%rng_stream_uniform%next()
    1300        1062 :                      rvek(ic) = 2.0_dp*rvek(ic) - 1.0_dp
    1301        1416 :                      rvek(ic) = rvek(ic)*helium_env(k)%helium%droplet_radius
    1302             :                   END DO
    1303        1416 :                   centroids(:, ia) = rvek(:)
    1304             :                   ! check if helium is outside of the droplet
    1305        1416 :                   r2 = DOT_PRODUCT(rvek, rvek)
    1306         354 :                   IF (r2 > helium_env(k)%helium%droplet_radius**2) THEN
    1307             :                      invalidpos = .TRUE.
    1308             :                   ELSE
    1309             :                      ! check for helium-helium collision
    1310        2340 :                      DO id = 1, ia - 1
    1311        8648 :                         rvek = centroids(:, ia) - centroids(:, id)
    1312        8648 :                         r2 = DOT_PRODUCT(rvek, rvek)
    1313        2340 :                         IF (r2 < minHeHedsttmp**2) THEN
    1314             :                            invalidpos = .TRUE.
    1315             :                            EXIT
    1316             :                         END IF
    1317             :                      END DO
    1318             :                   END IF
    1319         364 :                   IF (.NOT. invalidpos) THEN
    1320             :                      ! make sure the helium does not collide with the solute
    1321         178 :                      IF (helium_env(k)%helium%solute_present) THEN
    1322         712 :                         rvek(:) = (cvek(:) - centroids(:, ia))
    1323         712 :                         r2 = DOT_PRODUCT(rvek, rvek)
    1324         178 :                         IF (r2 <= solute_radius**2) invalidpos = .TRUE.
    1325             :                      END IF
    1326             :                   END IF
    1327             :                END DO
    1328             :             END DO
    1329             :             ! do thermal gaussian delocalization of hot start
    1330          14 :             IF (initkT > 0.0_dp) THEN
    1331           5 :                CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkT)
    1332             :             ELSE
    1333         162 :                DO ia = 1, helium_env(k)%helium%atoms
    1334        2610 :                   DO ib = 1, helium_env(k)%helium%beads
    1335        9945 :                      helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia)
    1336             :                   END DO
    1337             :                END DO
    1338             :             END IF
    1339         192 :             DO ia = 1, helium_env(k)%helium%atoms
    1340        3040 :                DO ib = 1, helium_env(k)%helium%beads
    1341             :                   ! Make sure, that nothing lies outside the droplet radius
    1342             :                   r1 = DOT_PRODUCT(helium_env(k)%helium%pos(:, ia, ib), &
    1343       11392 :                                    helium_env(k)%helium%pos(:, ia, ib))
    1344        2848 :                   IF (r1 > helium_env(k)%helium%droplet_radius**2) THEN
    1345          16 :                      r1 = SQRT(r1)
    1346             :                      helium_env(k)%helium%pos(:, ia, ib) = &
    1347             :                         helium_env(k)%helium%droplet_radius/r1* &
    1348          64 :                         helium_env(k)%helium%pos(:, ia, ib)
    1349        2832 :                   ELSE IF (helium_env(k)%helium%solute_present) THEN
    1350        2832 :                      IF (r1 < solute_radius**2) THEN
    1351             :                         !make sure that nothing lies within the molecule
    1352          16 :                         r1 = SQRT(r1)
    1353             :                         helium_env(k)%helium%pos(:, ia, ib) = &
    1354             :                            solute_radius/r1* &
    1355          64 :                            helium_env(k)%helium%pos(:, ia, ib)
    1356             :                      END IF
    1357             :                   END IF
    1358             :                   ! transfer to position around actual center of droplet
    1359             :                   helium_env(k)%helium%pos(:, ia, ib) = &
    1360             :                      helium_env(k)%helium%pos(:, ia, ib) + &
    1361       11570 :                      helium_env(k)%helium%center(:)
    1362             :                END DO
    1363             :             END DO
    1364             :          END IF
    1365       39246 :          helium_env(k)%helium%work = helium_env(k)%helium%pos
    1366          42 :          DEALLOCATE (centroids)
    1367             :       END DO
    1368             : 
    1369          18 :       RETURN
    1370          18 :    END SUBROUTINE helium_coord_init
    1371             : 
    1372             : ! **************************************************************************************************
    1373             : !> \brief ...
    1374             : !> \param helium_env ...
    1375             : !> \param centroids ...
    1376             : !> \param kbT ...
    1377             : ! **************************************************************************************************
    1378          15 :    SUBROUTINE helium_thermal_gaussian_beads_init(helium_env, centroids, kbT)
    1379             : 
    1380             :       TYPE(helium_solvent_type), POINTER                 :: helium_env
    1381             :       REAL(KIND=dp), DIMENSION(3, helium_env%atoms), &
    1382             :          INTENT(IN)                                      :: centroids
    1383             :       REAL(KIND=dp), INTENT(IN)                          :: kbT
    1384             : 
    1385             :       INTEGER                                            :: i, iatom, idim, imode, j, p
    1386             :       REAL(KIND=dp)                                      :: invsqrtp, omega, pip, rand, sqrt2p, &
    1387             :                                                             sqrtp, twopip, variance
    1388             :       REAL(KIND=dp), &
    1389          30 :          DIMENSION(helium_env%beads, helium_env%beads)   :: u2x
    1390          30 :       REAL(KIND=dp), DIMENSION(helium_env%beads)         :: nmhecoords
    1391             : 
    1392          15 :       p = helium_env%beads
    1393             : 
    1394          15 :       sqrt2p = SQRT(2.0_dp/REAL(p, dp))
    1395          15 :       twopip = twopi/REAL(p, dp)
    1396          15 :       pip = pi/REAL(p, dp)
    1397          15 :       sqrtp = SQRT(REAL(p, dp))
    1398          15 :       invsqrtp = 1.0_dp/SQRT(REAL(p, dp))
    1399             : 
    1400             :       ! set up normal mode backtransform matrix
    1401        6363 :       u2x(:, :) = 0.0_dp
    1402         309 :       u2x(:, 1) = invsqrtp
    1403         159 :       DO i = 2, p/2 + 1
    1404        3111 :          DO j = 1, p
    1405        3096 :             u2x(j, i) = sqrt2p*COS(twopip*(i - 1)*(j - 1))
    1406             :          END DO
    1407             :       END DO
    1408         150 :       DO i = p/2 + 2, p
    1409        2958 :          DO j = 1, p
    1410        2943 :             u2x(j, i) = sqrt2p*SIN(twopip*(i - 1)*(j - 1))
    1411             :          END DO
    1412             :       END DO
    1413          15 :       IF (MOD(p, 2) == 0) THEN
    1414           9 :          DO i = 1, p - 1, 2
    1415          72 :             u2x(i, p/2 + 1) = invsqrtp
    1416          72 :             u2x(i + 1, p/2 + 1) = -1.0_dp*invsqrtp
    1417             :          END DO
    1418             :       END IF
    1419             : 
    1420         360 :       DO iatom = 1, helium_env%atoms
    1421        1395 :          DO idim = 1, 3
    1422        1035 :             nmhecoords(1) = sqrtp*centroids(idim, iatom)
    1423       21744 :             DO imode = 2, p
    1424       20709 :                omega = 2.0_dp*p*kbT*SIN((imode - 1)*pip)
    1425       20709 :                variance = kbT*p/(helium_env%he_mass_au*omega**2)
    1426       20709 :                rand = helium_env%rng_stream_gaussian%next()
    1427       21744 :                nmhecoords(imode) = rand*SQRT(variance)
    1428             :             END DO
    1429      522372 :             helium_env%pos(idim, iatom, 1:p) = MATMUL(u2x, nmhecoords)
    1430             :          END DO
    1431             :       END DO
    1432             : 
    1433          15 :    END SUBROUTINE helium_thermal_gaussian_beads_init
    1434             : 
    1435             : ! ***************************************************************************
    1436             : !> \brief  Restore coordinates from the input structure.
    1437             : !> \param helium_env ...
    1438             : !> \date   2009-11-09
    1439             : !> \par    History
    1440             : !>         2010-07-22 accommodate additional cpus in the runtime wrt the
    1441             : !>                    restart [lwalewski]
    1442             : !>         2016-07-14 Modified to work with independent helium_env
    1443             : !>                    [cschran]
    1444             : !> \author Lukasz Walewski
    1445             : ! **************************************************************************************************
    1446           6 :    SUBROUTINE helium_coord_restore(helium_env)
    1447             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1448             : 
    1449             :       CHARACTER(len=default_string_length)               :: err_str, stmp
    1450             :       INTEGER                                            :: actlen, i, k, msglen, num_env_restart, &
    1451             :                                                             off, offset
    1452           6 :       LOGICAL, DIMENSION(:, :, :), POINTER               :: m
    1453           6 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: message
    1454           6 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: f
    1455             :       TYPE(cp_logger_type), POINTER                      :: logger
    1456             : 
    1457           6 :       NULLIFY (logger)
    1458          12 :       logger => cp_get_default_logger()
    1459             : 
    1460             :       ! assign the pointer to the memory location of the input structure, where
    1461             :       ! the coordinates are stored
    1462           6 :       NULLIFY (message)
    1463             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1464             :                                 "MOTION%PINT%HELIUM%COORD%_DEFAULT_KEYWORD_", &
    1465           6 :                                 r_vals=message)
    1466             : 
    1467             :       ! check that the number of values in the input match the current runtime
    1468           6 :       actlen = SIZE(message)
    1469           6 :       num_env_restart = actlen/helium_env(1)%helium%atoms/helium_env(1)%helium%beads/3
    1470             : 
    1471           6 :       IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
    1472           0 :          err_str = "Reading bead coordinates from the input file."
    1473           0 :          CALL helium_write_line(err_str)
    1474           0 :          err_str = "Number of environments in the restart...: '"
    1475           0 :          stmp = ""
    1476           0 :          WRITE (stmp, *) num_env_restart
    1477           0 :          err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
    1478           0 :          CALL helium_write_line(err_str)
    1479           0 :          err_str = "Number of current run time environments.: '"
    1480           0 :          stmp = ""
    1481           0 :          WRITE (stmp, *) helium_env(1)%helium%num_env
    1482           0 :          err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
    1483           0 :          CALL helium_write_line(err_str)
    1484           0 :          err_str = "Missmatch between number of bead coord. in input file and helium environments."
    1485           0 :          CPABORT(err_str)
    1486             :       ELSE
    1487           6 :          CALL helium_write_line("Bead coordinates read from the input file.")
    1488             : 
    1489           6 :          offset = 0
    1490           9 :          DO i = 1, logger%para_env%mepos
    1491           9 :             offset = offset + helium_env(1)%env_all(i)
    1492             :          END DO
    1493             : 
    1494             :          ! distribute coordinates over processors (no message passing)
    1495          12 :          DO k = 1, SIZE(helium_env)
    1496           6 :             msglen = helium_env(k)%helium%atoms*helium_env(k)%helium%beads*3
    1497           6 :             off = msglen*MOD(offset + k - 1, num_env_restart)
    1498             :             NULLIFY (m, f)
    1499          24 :             ALLOCATE (m(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads))
    1500          24 :             ALLOCATE (f(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads))
    1501       17034 :             m(:, :, :) = .TRUE.
    1502       17034 :             f(:, :, :) = 0.0_dp
    1503       17034 :             helium_env(k)%helium%pos(:, :, 1:helium_env(k)%helium%beads) = UNPACK(message(off + 1:off + msglen), MASK=m, FIELD=f)
    1504          12 :             DEALLOCATE (f, m)
    1505             :          END DO
    1506             : 
    1507             :       END IF
    1508             : 
    1509             :       NULLIFY (message)
    1510             : 
    1511           6 :       RETURN
    1512           6 :    END SUBROUTINE helium_coord_restore
    1513             : 
    1514             : ! ***************************************************************************
    1515             : !> \brief  Initialize forces exerted on the solute
    1516             : !> \param helium_env ...
    1517             : !> \date   2009-11-10
    1518             : !> \par    History
    1519             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1520             : !> \author Lukasz Walewski
    1521             : ! **************************************************************************************************
    1522           6 :    SUBROUTINE helium_force_init(helium_env)
    1523             : 
    1524             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1525             : 
    1526             :       INTEGER                                            :: k
    1527             : 
    1528          12 :       DO k = 1, SIZE(helium_env)
    1529          12 :          IF (helium_env(k)%helium%solute_present) THEN
    1530         780 :             helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
    1531         780 :             helium_env(k)%helium%force_inst(:, :) = 0.0_dp
    1532             :          END IF
    1533             :       END DO
    1534             : 
    1535           6 :       RETURN
    1536             :    END SUBROUTINE helium_force_init
    1537             : 
    1538             : ! ***************************************************************************
    1539             : !> \brief  Restore forces from the input structure to the runtime environment.
    1540             : !> \param helium_env ...
    1541             : !> \date   2009-11-10
    1542             : !> \par    History
    1543             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1544             : !> \author Lukasz Walewski
    1545             : ! **************************************************************************************************
    1546           2 :    SUBROUTINE helium_force_restore(helium_env)
    1547             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1548             : 
    1549             :       CHARACTER(len=default_string_length)               :: err_str, stmp
    1550             :       INTEGER                                            :: actlen, k, msglen
    1551           2 :       LOGICAL, DIMENSION(:, :), POINTER                  :: m
    1552           2 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: message
    1553           2 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: f
    1554             : 
    1555             : ! assign the pointer to the memory location of the input structure, where
    1556             : ! the forces are stored
    1557             : 
    1558           2 :       NULLIFY (message)
    1559             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1560             :                                 "MOTION%PINT%HELIUM%FORCE%_DEFAULT_KEYWORD_", &
    1561           2 :                                 r_vals=message)
    1562             : 
    1563             :       ! check if the destination array has correct size
    1564           2 :       msglen = helium_env(1)%helium%solute_atoms*helium_env(1)%helium%solute_beads*3
    1565           2 :       actlen = SIZE(helium_env(1)%helium%force_avrg)
    1566           2 :       err_str = "Invalid size of helium%force_avrg array: actual '"
    1567           2 :       stmp = ""
    1568           2 :       WRITE (stmp, *) actlen
    1569             :       err_str = TRIM(ADJUSTL(err_str))// &
    1570           2 :                 TRIM(ADJUSTL(stmp))//"' but expected '"
    1571           2 :       stmp = ""
    1572           2 :       WRITE (stmp, *) msglen
    1573           2 :       IF (actlen /= msglen) THEN
    1574             :          err_str = TRIM(ADJUSTL(err_str))// &
    1575           0 :                    TRIM(ADJUSTL(stmp))//"'."
    1576           0 :          CPABORT(err_str)
    1577             :       END IF
    1578             : 
    1579             :       ! restore forces on all processors (no message passing)
    1580             :       NULLIFY (m, f)
    1581           8 :       ALLOCATE (m(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3))
    1582           8 :       ALLOCATE (f(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3))
    1583          92 :       m(:, :) = .TRUE.
    1584          92 :       f(:, :) = 0.0_dp
    1585           4 :       DO k = 1, SIZE(helium_env)
    1586          92 :          helium_env(k)%helium%force_avrg(:, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f)
    1587          94 :          helium_env(k)%helium%force_inst(:, :) = 0.0_dp
    1588             :       END DO
    1589           2 :       DEALLOCATE (f, m)
    1590             : 
    1591           2 :       CALL helium_write_line("Forces on the solute read from the input file.")
    1592             : 
    1593             :       NULLIFY (message)
    1594             : 
    1595           2 :       RETURN
    1596           2 :    END SUBROUTINE helium_force_restore
    1597             : 
    1598             : ! ***************************************************************************
    1599             : !> \brief  Initialize the permutation state.
    1600             : !> \param helium_env ...
    1601             : !> \date   2009-11-05
    1602             : !> \par    History
    1603             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1604             : !> \author Lukasz Walewski
    1605             : !> \note   Assign the identity permutation at each processor. Inverse
    1606             : !>         permutation array gets assigned as well.
    1607             : ! **************************************************************************************************
    1608          18 :    SUBROUTINE helium_perm_init(helium_env)
    1609             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1610             : 
    1611             :       INTEGER                                            :: ia, k
    1612             : 
    1613          42 :       DO k = 1, SIZE(helium_env)
    1614         540 :          DO ia = 1, helium_env(k)%helium%atoms
    1615         498 :             helium_env(k)%helium%permutation(ia) = ia
    1616         522 :             helium_env(k)%helium%iperm(ia) = ia
    1617             :          END DO
    1618             :       END DO
    1619             : 
    1620          18 :       RETURN
    1621             :    END SUBROUTINE helium_perm_init
    1622             : 
    1623             : ! ***************************************************************************
    1624             : !> \brief  Restore permutation state from the input structure.
    1625             : !> \param helium_env ...
    1626             : !> \date   2009-11-05
    1627             : !> \par    History
    1628             : !>         2010-07-22 accommodate additional cpus in the runtime wrt the
    1629             : !>                    restart [lwalewski]
    1630             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1631             : !> \author Lukasz Walewski
    1632             : !> \note   Transfer permutation state from the input tree to the runtime
    1633             : !>         data structures on each processor. Inverse permutation array is
    1634             : !>         recalculated according to the restored permutation state.
    1635             : ! **************************************************************************************************
    1636           6 :    SUBROUTINE helium_perm_restore(helium_env)
    1637             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1638             : 
    1639             :       CHARACTER(len=default_string_length)               :: err_str, stmp
    1640             :       INTEGER                                            :: actlen, i, ia, ic, k, msglen, &
    1641             :                                                             num_env_restart, off, offset
    1642           6 :       INTEGER, DIMENSION(:), POINTER                     :: message
    1643             :       TYPE(cp_logger_type), POINTER                      :: logger
    1644             : 
    1645           6 :       NULLIFY (logger)
    1646          12 :       logger => cp_get_default_logger()
    1647             : 
    1648             :       ! assign the pointer to the memory location of the input structure, where
    1649             :       ! the permutation state is stored
    1650           6 :       NULLIFY (message)
    1651             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1652             :                                 "MOTION%PINT%HELIUM%PERM%_DEFAULT_KEYWORD_", &
    1653           6 :                                 i_vals=message)
    1654             : 
    1655             :       ! check the number of environments presumably stored in the restart
    1656           6 :       actlen = SIZE(message)
    1657           6 :       num_env_restart = actlen/helium_env(1)%helium%atoms
    1658             : 
    1659             : !TODO maybe add some sanity checks here:
    1660             : ! is num_env_restart integer ?
    1661             : 
    1662           6 :       IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
    1663           0 :          err_str = "Reading permutation state from the input file."
    1664           0 :          CALL helium_write_line(err_str)
    1665           0 :          err_str = "Number of environments in the restart...: '"
    1666           0 :          stmp = ""
    1667           0 :          WRITE (stmp, *) num_env_restart
    1668           0 :          err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
    1669           0 :          CALL helium_write_line(err_str)
    1670           0 :          err_str = "Number of current run time environments.: '"
    1671           0 :          stmp = ""
    1672           0 :          WRITE (stmp, *) helium_env(1)%helium%num_env
    1673           0 :          err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
    1674           0 :          CALL helium_write_line(err_str)
    1675           0 :          err_str = "Missmatch between number of perm. states in input file and helium environments."
    1676           0 :          CPABORT(err_str)
    1677             :       ELSE
    1678           6 :          CALL helium_write_line("Permutation state read from the input file.")
    1679             : 
    1680             :          ! distribute permutation state over processors
    1681           6 :          offset = 0
    1682           9 :          DO i = 1, logger%para_env%mepos
    1683           9 :             offset = offset + helium_env(1)%env_all(i)
    1684             :          END DO
    1685             : 
    1686          12 :          DO k = 1, SIZE(helium_env)
    1687           6 :             msglen = helium_env(k)%helium%atoms
    1688           6 :             off = msglen*MOD(k - 1 + offset, num_env_restart)
    1689         402 :             helium_env(k)%helium%permutation(:) = message(off + 1:off + msglen)
    1690             :          END DO
    1691             :       END IF
    1692             : 
    1693             :       ! recalculate the inverse permutation array
    1694          12 :       DO k = 1, SIZE(helium_env)
    1695         198 :          helium_env(k)%helium%iperm(:) = 0
    1696           6 :          ic = 0
    1697         198 :          DO ia = 1, msglen
    1698         198 :             IF ((helium_env(k)%helium%permutation(ia) > 0) .AND. (helium_env(k)%helium%permutation(ia) <= msglen)) THEN
    1699         192 :                helium_env(k)%helium%iperm(helium_env(k)%helium%permutation(ia)) = ia
    1700         192 :                ic = ic + 1
    1701             :             END IF
    1702             :          END DO
    1703           6 :          err_str = "Invalid HELIUM%PERM state: some numbers not within (1,"
    1704           6 :          stmp = ""
    1705           6 :          WRITE (stmp, *) msglen
    1706          12 :          IF (ic /= msglen) THEN
    1707             :             err_str = TRIM(ADJUSTL(err_str))// &
    1708           0 :                       TRIM(ADJUSTL(stmp))//")."
    1709           0 :             CPABORT(err_str)
    1710             :          END IF
    1711             :       END DO
    1712             :       NULLIFY (message)
    1713             : 
    1714           6 :       RETURN
    1715           6 :    END SUBROUTINE helium_perm_restore
    1716             : 
    1717             : ! ***************************************************************************
    1718             : !> \brief  Restore averages from the input structure
    1719             : !> \param helium_env ...
    1720             : !> \date   2014-06-25
    1721             : !> \par    History
    1722             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1723             : !> \author Lukasz Walewski
    1724             : ! **************************************************************************************************
    1725         120 :    SUBROUTINE helium_averages_restore(helium_env)
    1726             : 
    1727             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1728             : 
    1729             :       INTEGER                                            :: i, k, msglen, num_env_restart, off, &
    1730             :                                                             offset
    1731             :       LOGICAL                                            :: explicit
    1732          24 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: message
    1733             :       TYPE(cp_logger_type), POINTER                      :: logger
    1734             : 
    1735          24 :       NULLIFY (logger)
    1736          48 :       logger => cp_get_default_logger()
    1737             : 
    1738          24 :       offset = 0
    1739          36 :       DO i = 1, logger%para_env%mepos
    1740          36 :          offset = offset + helium_env(1)%env_all(i)
    1741             :       END DO
    1742             : 
    1743             :       ! restore projected area
    1744             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1745             :                                 "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
    1746          24 :                                 explicit=explicit)
    1747          24 :       IF (explicit) THEN
    1748           0 :          NULLIFY (message)
    1749             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    1750             :                                    "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
    1751           0 :                                    r_vals=message)
    1752           0 :          num_env_restart = SIZE(message)/3 ! apparent number of environments
    1753           0 :          msglen = 3
    1754           0 :          DO k = 1, SIZE(helium_env)
    1755           0 :             off = msglen*MOD(offset + k - 1, num_env_restart)
    1756           0 :             helium_env(k)%helium%proarea%rstr(:) = message(off + 1:off + msglen)
    1757             :          END DO
    1758             :       ELSE
    1759          54 :          DO k = 1, SIZE(helium_env)
    1760         144 :             helium_env(k)%helium%proarea%rstr(:) = 0.0_dp
    1761             :          END DO
    1762             :       END IF
    1763             : 
    1764             :       ! restore projected area squared
    1765             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1766             :                                 "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
    1767          24 :                                 explicit=explicit)
    1768          24 :       IF (explicit) THEN
    1769           0 :          NULLIFY (message)
    1770             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    1771             :                                    "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
    1772           0 :                                    r_vals=message)
    1773           0 :          num_env_restart = SIZE(message)/3 ! apparent number of environments
    1774           0 :          msglen = 3
    1775           0 :          DO k = 1, SIZE(helium_env)
    1776           0 :             off = msglen*MOD(offset + k - 1, num_env_restart)
    1777           0 :             helium_env(k)%helium%prarea2%rstr(:) = message(off + 1:off + msglen)
    1778             :          END DO
    1779             :       ELSE
    1780          54 :          DO k = 1, SIZE(helium_env)
    1781         144 :             helium_env(k)%helium%prarea2%rstr(:) = 0.0_dp
    1782             :          END DO
    1783             :       END IF
    1784             : 
    1785             :       ! restore winding number squared
    1786             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1787             :                                 "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
    1788          24 :                                 explicit=explicit)
    1789          24 :       IF (explicit) THEN
    1790           0 :          NULLIFY (message)
    1791             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    1792             :                                    "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
    1793           0 :                                    r_vals=message)
    1794           0 :          num_env_restart = SIZE(message)/3 ! apparent number of environments
    1795           0 :          msglen = 3
    1796           0 :          DO k = 1, SIZE(helium_env)
    1797           0 :             off = msglen*MOD(offset + k - 1, num_env_restart)
    1798           0 :             helium_env(k)%helium%wnmber2%rstr(:) = message(off + 1:off + msglen)
    1799             :          END DO
    1800             :       ELSE
    1801          54 :          DO k = 1, SIZE(helium_env)
    1802         144 :             helium_env(k)%helium%wnmber2%rstr(:) = 0.0_dp
    1803             :          END DO
    1804             :       END IF
    1805             : 
    1806             :       ! restore moment of inertia
    1807             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1808             :                                 "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
    1809          24 :                                 explicit=explicit)
    1810          24 :       IF (explicit) THEN
    1811           0 :          NULLIFY (message)
    1812             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    1813             :                                    "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
    1814           0 :                                    r_vals=message)
    1815           0 :          num_env_restart = SIZE(message)/3 ! apparent number of environments
    1816           0 :          msglen = 3
    1817           0 :          DO k = 1, SIZE(helium_env)
    1818           0 :             off = msglen*MOD(offset + k - 1, num_env_restart)
    1819           0 :             helium_env(k)%helium%mominer%rstr(:) = message(off + 1:off + msglen)
    1820             :          END DO
    1821             :       ELSE
    1822          54 :          DO k = 1, SIZE(helium_env)
    1823         144 :             helium_env(k)%helium%mominer%rstr(:) = 0.0_dp
    1824             :          END DO
    1825             :       END IF
    1826             : 
    1827          24 :       IF (helium_env(1)%helium%rdf_present) THEN
    1828          16 :          CALL helium_rdf_restore(helium_env)
    1829             :       END IF
    1830             : 
    1831          24 :       IF (helium_env(1)%helium%rho_present) THEN
    1832             :          ! restore densities
    1833          24 :          CALL helium_rho_restore(helium_env)
    1834             :       END IF
    1835             : 
    1836             :       ! get the weighting factor
    1837          54 :       DO k = 1, SIZE(helium_env)
    1838             :          CALL section_vals_val_get( &
    1839             :             helium_env(k)%helium%input, &
    1840             :             "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
    1841          30 :             i_val=helium_env(k)%helium%averages_iweight)
    1842             : 
    1843             :          ! set the flag indicating whether the averages have been restarted
    1844             :          CALL section_vals_val_get( &
    1845             :             helium_env(k)%helium%input, &
    1846             :             "EXT_RESTART%RESTART_HELIUM_AVERAGES", &
    1847          54 :             l_val=helium_env(k)%helium%averages_restarted)
    1848             :       END DO
    1849             : 
    1850          24 :       RETURN
    1851          24 :    END SUBROUTINE helium_averages_restore
    1852             : 
    1853             : ! ***************************************************************************
    1854             : !> \brief  Create RNG streams and initialize their state.
    1855             : !> \param helium_env ...
    1856             : !> \date   2009-11-04
    1857             : !> \par    History
    1858             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1859             : !> \author Lukasz Walewski
    1860             : !> \note   TODO: This function shouldn't create (allocate) objects! Only
    1861             : !>         initialization, i.e. setting the seed values etc, should be done
    1862             : !>         here, allocation should be moved to helium_create
    1863             : ! **************************************************************************************************
    1864          24 :    SUBROUTINE helium_rng_init(helium_env)
    1865             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1866             : 
    1867             :       INTEGER                                            :: helium_seed, i, offset
    1868             :       REAL(KIND=dp), DIMENSION(3, 2)                     :: initial_seed
    1869             :       TYPE(cp_logger_type), POINTER                      :: logger
    1870          24 :       TYPE(rng_stream_p_type), DIMENSION(:), POINTER     :: gaussian_array, uniform_array
    1871             : 
    1872          24 :       NULLIFY (logger)
    1873          48 :       logger => cp_get_default_logger()
    1874          24 :       IF (logger%para_env%is_source()) THEN
    1875             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    1876             :                                    "MOTION%PINT%HELIUM%RNG_SEED", &
    1877          12 :                                    i_val=helium_seed)
    1878             :       END IF
    1879             :       CALL helium_env(1)%comm%bcast(helium_seed, &
    1880          24 :                                     logger%para_env%source)
    1881         216 :       initial_seed(:, :) = REAL(helium_seed, dp)
    1882             : 
    1883             :       ALLOCATE (uniform_array(helium_env(1)%helium%num_env), &
    1884         240 :                 gaussian_array(helium_env(1)%helium%num_env))
    1885          84 :       DO i = 1, helium_env(1)%helium%num_env
    1886        2964 :          ALLOCATE (uniform_array(i)%stream, gaussian_array(i)%stream)
    1887             :       END DO
    1888             : 
    1889             :       ! Create num_env RNG streams on processor all processors
    1890             :       ! and distribute them so, that each processor gets unique
    1891             :       ! RN sequences for his helium environments
    1892             :       ! COMMENT: rng_stream can not be used with mp_bcast
    1893             : 
    1894             :       uniform_array(1)%stream = rng_stream_type(name="helium_rns_uniform", &
    1895             :                                                 distribution_type=UNIFORM, &
    1896             :                                                 extended_precision=.TRUE., &
    1897          24 :                                                 seed=initial_seed)
    1898             : 
    1899             :       gaussian_array(1)%stream = rng_stream_type(name="helium_rns_gaussian", &
    1900             :                                                  distribution_type=GAUSSIAN, &
    1901             :                                                  extended_precision=.TRUE., &
    1902          24 :                                                  last_rng_stream=uniform_array(1)%stream)
    1903          60 :       DO i = 2, helium_env(1)%helium%num_env
    1904             :          uniform_array(i)%stream = rng_stream_type(name="helium_rns_uniform", &
    1905             :                                                    distribution_type=UNIFORM, &
    1906             :                                                    extended_precision=.TRUE., &
    1907          36 :                                                    last_rng_stream=gaussian_array(i - 1)%stream)
    1908             : 
    1909             :          gaussian_array(i)%stream = rng_stream_type(name="helium_rns_uniform", &
    1910             :                                                     distribution_type=GAUSSIAN, &
    1911             :                                                     extended_precision=.TRUE., &
    1912          60 :                                                     last_rng_stream=uniform_array(i)%stream)
    1913             :       END DO
    1914             : 
    1915          24 :       offset = 0
    1916          36 :       DO i = 1, logger%para_env%mepos
    1917          36 :          offset = offset + helium_env(1)%env_all(i)
    1918             :       END DO
    1919             : 
    1920          54 :       DO i = 1, SIZE(helium_env)
    1921             :          NULLIFY (helium_env(i)%helium%rng_stream_uniform, &
    1922          30 :                   helium_env(i)%helium%rng_stream_gaussian)
    1923          30 :          helium_env(i)%helium%rng_stream_uniform => uniform_array(offset + i)%stream
    1924          54 :          helium_env(i)%helium%rng_stream_gaussian => gaussian_array(offset + i)%stream
    1925             :       END DO
    1926             : 
    1927          84 :       DO i = 1, helium_env(1)%helium%num_env
    1928          60 :          IF (i .LE. offset .OR. i .GT. offset + SIZE(helium_env)) THEN
    1929             :             ! only deallocate pointers here which were not passed on to helium_env(*)%helium
    1930          30 :             DEALLOCATE (uniform_array(i)%stream)
    1931          30 :             DEALLOCATE (gaussian_array(i)%stream)
    1932             :          END IF
    1933          60 :          NULLIFY (uniform_array(i)%stream)
    1934          84 :          NULLIFY (gaussian_array(i)%stream)
    1935             :       END DO
    1936             : 
    1937          24 :       DEALLOCATE (uniform_array)
    1938          24 :       DEALLOCATE (gaussian_array)
    1939          24 :    END SUBROUTINE helium_rng_init
    1940             : 
    1941             : ! ***************************************************************************
    1942             : !> \brief  Restore RNG state from the input structure.
    1943             : !> \param helium_env ...
    1944             : !> \date   2009-11-04
    1945             : !> \par    History
    1946             : !>         2010-07-22 Create new rng streams if more cpus available in the
    1947             : !>         runtime than in the restart [lwalewski]
    1948             : !>         2016-04-18 Modified for independet number of helium_env [cschran]
    1949             : !> \author Lukasz Walewski
    1950             : ! **************************************************************************************************
    1951           6 :    SUBROUTINE helium_rng_restore(helium_env)
    1952             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1953             : 
    1954             :       CHARACTER(len=default_string_length)               :: err_str, stmp
    1955             :       INTEGER                                            :: actlen, i, k, msglen, num_env_restart, &
    1956             :                                                             off, offset
    1957             :       LOGICAL                                            :: lbf
    1958             :       LOGICAL, DIMENSION(3, 2)                           :: m
    1959             :       REAL(KIND=dp)                                      :: bf, bu
    1960             :       REAL(KIND=dp), DIMENSION(3, 2)                     :: bg, cg, f, ig
    1961           6 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: message
    1962             :       TYPE(cp_logger_type), POINTER                      :: logger
    1963             : 
    1964           6 :       NULLIFY (logger)
    1965          12 :       logger => cp_get_default_logger()
    1966             : 
    1967             :       ! assign the pointer to the memory location of the input structure
    1968             :       ! where the RNG state is stored
    1969           6 :       NULLIFY (message)
    1970             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    1971             :                                 "MOTION%PINT%HELIUM%RNG_STATE%_DEFAULT_KEYWORD_", &
    1972           6 :                                 r_vals=message)
    1973             : 
    1974             :       ! check the number of environments presumably stored in the restart
    1975           6 :       actlen = SIZE(message)
    1976           6 :       num_env_restart = actlen/40
    1977             : 
    1978             :       ! check, if RNG restart has the same dimension as helium%num_env
    1979           6 :       IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
    1980           0 :          err_str = "Reading RNG state from the input file."
    1981           0 :          CALL helium_write_line(err_str)
    1982           0 :          err_str = "Number of environments in the restart...: '"
    1983           0 :          stmp = ""
    1984           0 :          WRITE (stmp, *) num_env_restart
    1985           0 :          err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
    1986           0 :          CALL helium_write_line(err_str)
    1987           0 :          err_str = "Number of current run time environments.: '"
    1988           0 :          stmp = ""
    1989           0 :          WRITE (stmp, *) helium_env(1)%helium%num_env
    1990           0 :          err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
    1991           0 :          CALL helium_write_line(err_str)
    1992           0 :          err_str = "Missmatch between number of RNG states in input file and helium environments."
    1993           0 :          CPABORT(err_str)
    1994             :       ELSE
    1995           6 :          CALL helium_write_line("RNG state read from the input file.")
    1996             : 
    1997             :          ! unpack the buffer at each processor, set RNG state
    1998           6 :          offset = 0
    1999           9 :          DO i = 1, logger%para_env%mepos
    2000           9 :             offset = offset + helium_env(1)%env_all(i)
    2001             :          END DO
    2002             : 
    2003          12 :          DO k = 1, SIZE(helium_env)
    2004           6 :             msglen = 40
    2005           6 :             off = msglen*(offset + k - 1)
    2006          54 :             m(:, :) = .TRUE.
    2007           6 :             f(:, :) = 0.0_dp
    2008           6 :             bg(:, :) = UNPACK(message(off + 1:off + 6), MASK=m, FIELD=f)
    2009           6 :             cg(:, :) = UNPACK(message(off + 7:off + 12), MASK=m, FIELD=f)
    2010           6 :             ig(:, :) = UNPACK(message(off + 13:off + 18), MASK=m, FIELD=f)
    2011           6 :             bf = message(off + 19)
    2012           6 :             bu = message(off + 20)
    2013           6 :             IF (bf .GT. 0) THEN
    2014           0 :                lbf = .TRUE.
    2015             :             ELSE
    2016           6 :                lbf = .FALSE.
    2017             :             END IF
    2018             :             CALL helium_env(k)%helium%rng_stream_uniform%set(bg=bg, cg=cg, ig=ig, &
    2019           6 :                                                              buffer=bu, buffer_filled=lbf)
    2020           6 :             bg(:, :) = UNPACK(message(off + 21:off + 26), MASK=m, FIELD=f)
    2021           6 :             cg(:, :) = UNPACK(message(off + 27:off + 32), MASK=m, FIELD=f)
    2022           6 :             ig(:, :) = UNPACK(message(off + 33:off + 38), MASK=m, FIELD=f)
    2023           6 :             bf = message(off + 39)
    2024           6 :             bu = message(off + 40)
    2025           6 :             IF (bf .GT. 0) THEN
    2026           4 :                lbf = .TRUE.
    2027             :             ELSE
    2028           2 :                lbf = .FALSE.
    2029             :             END IF
    2030             :             CALL helium_env(k)%helium%rng_stream_gaussian%set(bg=bg, cg=cg, ig=ig, &
    2031          12 :                                                               buffer=bu, buffer_filled=lbf)
    2032             :          END DO
    2033             :       END IF
    2034             : 
    2035             :       NULLIFY (message)
    2036             : 
    2037           6 :       RETURN
    2038           6 :    END SUBROUTINE helium_rng_restore
    2039             : 
    2040             : ! ***************************************************************************
    2041             : !> \brief  Create the RDF related data structures and initialize
    2042             : !> \param helium ...
    2043             : !> \date   2014-02-25
    2044             : !> \par    History
    2045             : !>         2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
    2046             : !> \author Lukasz Walewski
    2047             : ! **************************************************************************************************
    2048          30 :    SUBROUTINE helium_rdf_init(helium)
    2049             : 
    2050             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2051             : 
    2052             :       CHARACTER(len=2*default_string_length)             :: err_str, stmp
    2053             :       INTEGER                                            :: ii, ij
    2054             :       LOGICAL                                            :: explicit
    2055             :       TYPE(cp_logger_type), POINTER                      :: logger
    2056             : 
    2057             :       ! read parameters
    2058          30 :       NULLIFY (logger)
    2059          30 :       logger => cp_get_default_logger()
    2060             :       CALL section_vals_val_get( &
    2061             :          helium%input, &
    2062             :          "MOTION%PINT%HELIUM%RDF%SOLUTE_HE", &
    2063          30 :          l_val=helium%rdf_sol_he)
    2064             :       CALL section_vals_val_get( &
    2065             :          helium%input, &
    2066             :          "MOTION%PINT%HELIUM%RDF%HE_HE", &
    2067          30 :          l_val=helium%rdf_he_he)
    2068             : 
    2069             :       ! Prevent He_He Rdfs for a single he atom:
    2070          30 :       IF (helium%atoms <= 1) THEN
    2071           0 :          helium%rdf_he_he = .FALSE.
    2072             :       END IF
    2073             : 
    2074          30 :       helium%rdf_num = 0
    2075          30 :       IF (helium%rdf_sol_he) THEN
    2076          30 :          IF (helium%solute_present) THEN
    2077             :             ! get number of centers from solute to store solute positions
    2078          80 :             ALLOCATE (helium%rdf_centers(helium%beads, helium%solute_atoms*3))
    2079        3080 :             helium%rdf_centers(:, :) = 0.0_dp
    2080          20 :             helium%rdf_num = helium%solute_atoms
    2081             :          ELSE
    2082          10 :             helium%rdf_sol_he = .FALSE.
    2083             :          END IF
    2084             :       END IF
    2085             : 
    2086          30 :       IF (helium%rdf_he_he) THEN
    2087           2 :          helium%rdf_num = helium%rdf_num + 1
    2088             :       END IF
    2089             : 
    2090             :       ! set the flag for RDF and either proceed or return
    2091          30 :       IF (helium%rdf_num > 0) THEN
    2092          22 :          helium%rdf_present = .TRUE.
    2093             :       ELSE
    2094           8 :          helium%rdf_present = .FALSE.
    2095             :          err_str = "HELIUM RDFs requested, but no valid choice of "// &
    2096           8 :                    "parameters specified. No RDF will be computed."
    2097           8 :          CPWARN(err_str)
    2098           8 :          RETURN
    2099             :       END IF
    2100             : 
    2101             :       ! set the maximum RDF range
    2102             :       CALL section_vals_val_get( &
    2103             :          helium%input, &
    2104             :          "MOTION%PINT%HELIUM%RDF%MAXR", &
    2105          22 :          explicit=explicit)
    2106          22 :       IF (explicit) THEN
    2107             :          ! use the value explicitly set in the input
    2108             :          CALL section_vals_val_get( &
    2109             :             helium%input, &
    2110             :             "MOTION%PINT%HELIUM%RDF%MAXR", &
    2111           0 :             r_val=helium%rdf_maxr)
    2112             :       ELSE
    2113             :          ! use the default value
    2114             :          CALL section_vals_val_get( &
    2115             :             helium%input, &
    2116             :             "MOTION%PINT%HELIUM%DROPLET_RADIUS", &
    2117          22 :             explicit=explicit)
    2118          22 :          IF (explicit) THEN
    2119             :             ! use the droplet radius
    2120          14 :             IF (helium%solute_present .AND. .NOT. helium%periodic) THEN
    2121             :                ! COM of solute is used as center of the box.
    2122             :                ! Therefore distances became larger then droplet_radius
    2123             :                ! when parts of the solute are on opposite site of
    2124             :                ! COM and helium.
    2125             :                ! Use two times droplet_radius for security:
    2126          14 :                helium%rdf_maxr = helium%droplet_radius*2.0_dp
    2127             :             ELSE
    2128           0 :                helium%rdf_maxr = helium%droplet_radius
    2129             :             END IF
    2130             :          ELSE
    2131             :             ! use cell_size and cell_shape
    2132             :             ! (they are set regardless of us being periodic or not)
    2133           8 :             SELECT CASE (helium%cell_shape)
    2134             :             CASE (helium_cell_shape_cube)
    2135           0 :                helium%rdf_maxr = helium%cell_size/2.0_dp
    2136             :             CASE (helium_cell_shape_octahedron)
    2137           8 :                helium%rdf_maxr = helium%cell_size*SQRT(3.0_dp)/4.0_dp
    2138             :             CASE DEFAULT
    2139           0 :                helium%rdf_maxr = 0.0_dp
    2140           0 :                WRITE (stmp, *) helium%cell_shape
    2141           0 :                err_str = "cell shape unknown ("//TRIM(ADJUSTL(stmp))//")"
    2142           8 :                CPABORT(err_str)
    2143             :             END SELECT
    2144             :          END IF
    2145             :       END IF
    2146             : 
    2147             :       ! get number of bins and set bin spacing
    2148             :       CALL section_vals_val_get( &
    2149             :          helium%input, &
    2150             :          "MOTION%PINT%HELIUM%RDF%NBIN", &
    2151          22 :          i_val=helium%rdf_nbin)
    2152          22 :       helium%rdf_delr = helium%rdf_maxr/REAL(helium%rdf_nbin, dp)
    2153             : 
    2154             :       ! get the weighting factor
    2155             :       CALL section_vals_val_get( &
    2156             :          helium%input, &
    2157             :          "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
    2158          22 :          i_val=helium%rdf_iweight)
    2159             : 
    2160             :       ! allocate and initialize memory for RDF storage
    2161          22 :       ii = helium%rdf_num
    2162          22 :       ij = helium%rdf_nbin
    2163          88 :       ALLOCATE (helium%rdf_inst(ii, ij))
    2164          88 :       ALLOCATE (helium%rdf_accu(ii, ij))
    2165          88 :       ALLOCATE (helium%rdf_rstr(ii, ij))
    2166       21022 :       helium%rdf_inst(:, :) = 0.0_dp
    2167       21022 :       helium%rdf_accu(:, :) = 0.0_dp
    2168       21022 :       helium%rdf_rstr(:, :) = 0.0_dp
    2169             : 
    2170             :       RETURN
    2171          22 :    END SUBROUTINE helium_rdf_init
    2172             : 
    2173             : ! ***************************************************************************
    2174             : !> \brief  Restore the RDFs from the input structure
    2175             : !> \param helium_env ...
    2176             : !> \date   2011-06-22
    2177             : !> \par    History
    2178             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    2179             : !>         2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
    2180             : !> \author Lukasz Walewski
    2181             : ! **************************************************************************************************
    2182          16 :    SUBROUTINE helium_rdf_restore(helium_env)
    2183             : 
    2184             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    2185             : 
    2186             :       CHARACTER(len=default_string_length)               :: stmp1, stmp2
    2187             :       CHARACTER(len=max_line_length)                     :: err_str
    2188             :       INTEGER                                            :: ii, ij, itmp, k, msglen
    2189             :       LOGICAL                                            :: explicit, ltmp
    2190          16 :       LOGICAL, DIMENSION(:, :), POINTER                  :: m
    2191          16 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: message
    2192          16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: f
    2193             : 
    2194             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    2195             :                                 "MOTION%PINT%HELIUM%AVERAGES%RDF", &
    2196          16 :                                 explicit=explicit)
    2197          16 :       IF (explicit) THEN
    2198           0 :          NULLIFY (message)
    2199             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    2200             :                                    "MOTION%PINT%HELIUM%AVERAGES%RDF", &
    2201           0 :                                    r_vals=message)
    2202           0 :          msglen = SIZE(message)
    2203           0 :          itmp = SIZE(helium_env(1)%helium%rdf_rstr)
    2204           0 :          ltmp = (msglen .EQ. itmp)
    2205           0 :          IF (.NOT. ltmp) THEN
    2206           0 :             stmp1 = ""
    2207           0 :             WRITE (stmp1, *) msglen
    2208           0 :             stmp2 = ""
    2209           0 :             WRITE (stmp2, *) itmp
    2210             :             err_str = "Size of the RDF array in the input ("// &
    2211             :                       TRIM(ADJUSTL(stmp1))// &
    2212             :                       ") .NE. that in the runtime ("// &
    2213           0 :                       TRIM(ADJUSTL(stmp2))//")."
    2214           0 :             CPABORT(err_str)
    2215             :          END IF
    2216             :       ELSE
    2217             :          RETURN
    2218             :       END IF
    2219             : 
    2220           0 :       ii = helium_env(1)%helium%rdf_num
    2221           0 :       ij = helium_env(1)%helium%rdf_nbin
    2222             :       NULLIFY (m, f)
    2223           0 :       ALLOCATE (m(ii, ij))
    2224           0 :       ALLOCATE (f(ii, ij))
    2225           0 :       m(:, :) = .TRUE.
    2226           0 :       f(:, :) = 0.0_dp
    2227             : 
    2228           0 :       DO k = 1, SIZE(helium_env)
    2229           0 :          helium_env(k)%helium%rdf_rstr(:, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f)
    2230             :          CALL section_vals_val_get(helium_env(k)%helium%input, &
    2231             :                                    "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
    2232           0 :                                    i_val=helium_env(k)%helium%rdf_iweight)
    2233             :       END DO
    2234             : 
    2235           0 :       DEALLOCATE (f, m)
    2236             :       NULLIFY (message)
    2237             : 
    2238           0 :       RETURN
    2239          16 :    END SUBROUTINE helium_rdf_restore
    2240             : 
    2241             : ! ***************************************************************************
    2242             : !> \brief  Release/deallocate RDF related data structures
    2243             : !> \param helium ...
    2244             : !> \date   2014-02-25
    2245             : !> \author Lukasz Walewski
    2246             : ! **************************************************************************************************
    2247          30 :    SUBROUTINE helium_rdf_release(helium)
    2248             : 
    2249             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2250             : 
    2251          30 :       IF (helium%rdf_present) THEN
    2252             : 
    2253           0 :          DEALLOCATE ( &
    2254             :             helium%rdf_rstr, &
    2255           0 :             helium%rdf_accu, &
    2256          22 :             helium%rdf_inst)
    2257             : 
    2258             :          NULLIFY ( &
    2259             :             helium%rdf_rstr, &
    2260          22 :             helium%rdf_accu, &
    2261          22 :             helium%rdf_inst)
    2262             : 
    2263          22 :          IF (helium%rdf_sol_he) THEN
    2264          20 :             DEALLOCATE (helium%rdf_centers)
    2265          20 :             NULLIFY (helium%rdf_centers)
    2266             :          END IF
    2267             : 
    2268             :       END IF
    2269             : 
    2270          30 :       RETURN
    2271             :    END SUBROUTINE helium_rdf_release
    2272             : 
    2273             : ! ***************************************************************************
    2274             : !> \brief  Check whether property <prop> is activated in the input structure
    2275             : !> \param helium ...
    2276             : !> \param prop ...
    2277             : !> \return ...
    2278             : !> \date   2014-06-26
    2279             : !> \author Lukasz Walewski
    2280             : !> \note   The property is controlled by two items in the input structure,
    2281             : !>         the printkey and the control section. Two settings result in
    2282             : !>         the property being considered active:
    2283             : !>         1. printkey is on at the given print level
    2284             : !>         2. control section is explicit and on
    2285             : !>         If the property is considered active it should be calculated
    2286             : !>         and printed through out the run.
    2287             : ! **************************************************************************************************
    2288          64 :    FUNCTION helium_property_active(helium, prop) RESULT(is_active)
    2289             : 
    2290             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2291             :       CHARACTER(len=*)                                   :: prop
    2292             :       LOGICAL                                            :: is_active
    2293             : 
    2294             :       CHARACTER(len=default_string_length)               :: input_path
    2295             :       INTEGER                                            :: print_level
    2296             :       LOGICAL                                            :: explicit, is_on
    2297             :       TYPE(cp_logger_type), POINTER                      :: logger
    2298             :       TYPE(section_vals_type), POINTER                   :: print_key, section
    2299             : 
    2300          60 :       is_active = .FALSE.
    2301          60 :       NULLIFY (logger)
    2302         120 :       logger => cp_get_default_logger()
    2303             : 
    2304             :       ! if the printkey is on at this runlevel consider prop to be active
    2305          60 :       NULLIFY (print_key)
    2306          60 :       input_path = "MOTION%PINT%HELIUM%PRINT%"//TRIM(ADJUSTL(prop))
    2307             :       print_key => section_vals_get_subs_vals( &
    2308             :                    helium%input, &
    2309          60 :                    input_path)
    2310             :       is_on = cp_printkey_is_on( &
    2311             :               iteration_info=logger%iter_info, &
    2312          60 :               print_key=print_key)
    2313          60 :       IF (is_on) THEN
    2314          60 :          is_active = .TRUE.
    2315             :          RETURN
    2316             :       END IF
    2317             : 
    2318             :       ! if the control section is explicit and on consider prop to be active
    2319             :       ! and activate the printkey
    2320           4 :       is_active = .FALSE.
    2321           4 :       NULLIFY (section)
    2322           4 :       input_path = "MOTION%PINT%HELIUM%"//TRIM(ADJUSTL(prop))
    2323             :       section => section_vals_get_subs_vals( &
    2324             :                  helium%input, &
    2325           4 :                  input_path)
    2326           4 :       CALL section_vals_get(section, explicit=explicit)
    2327           4 :       IF (explicit) THEN
    2328             :          ! control section explicitly present, continue checking
    2329             :          CALL section_vals_val_get( &
    2330             :             section, &
    2331             :             "_SECTION_PARAMETERS_", &
    2332           4 :             l_val=is_on)
    2333           4 :          IF (is_on) THEN
    2334             :             ! control section is explicit and on, activate the property
    2335           4 :             is_active = .TRUE.
    2336             :             ! activate the corresponding print_level as well
    2337           4 :             print_level = logger%iter_info%print_level
    2338             :             CALL section_vals_val_set( &
    2339             :                print_key, &
    2340             :                "_SECTION_PARAMETERS_", &
    2341           4 :                i_val=print_level)
    2342             :          END IF
    2343             :       END IF
    2344             : 
    2345             :       RETURN
    2346          60 :    END FUNCTION helium_property_active
    2347             : 
    2348             : ! ***************************************************************************
    2349             : !> \brief  Create the density related data structures and initialize
    2350             : !> \param helium ...
    2351             : !> \date   2014-02-25
    2352             : !> \author Lukasz Walewski
    2353             : ! **************************************************************************************************
    2354          30 :    SUBROUTINE helium_rho_property_init(helium)
    2355             : 
    2356             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2357             : 
    2358             :       INTEGER                                            :: nc
    2359             : 
    2360          30 :       ALLOCATE (helium%rho_property(rho_num))
    2361             : 
    2362          30 :       helium%rho_property(rho_atom_number)%name = 'Atom number density'
    2363          30 :       nc = 1
    2364          30 :       helium%rho_property(rho_atom_number)%num_components = nc
    2365          30 :       ALLOCATE (helium%rho_property(rho_atom_number)%filename_suffix(nc))
    2366          30 :       ALLOCATE (helium%rho_property(rho_atom_number)%component_name(nc))
    2367          30 :       ALLOCATE (helium%rho_property(rho_atom_number)%component_index(nc))
    2368          30 :       helium%rho_property(rho_atom_number)%filename_suffix(1) = 'an'
    2369          30 :       helium%rho_property(rho_atom_number)%component_name(1) = ''
    2370          60 :       helium%rho_property(rho_atom_number)%component_index(:) = 0
    2371             : 
    2372          30 :       helium%rho_property(rho_projected_area)%name = 'Projected area squared density, A*A(r)'
    2373          30 :       nc = 3
    2374          30 :       helium%rho_property(rho_projected_area)%num_components = nc
    2375          30 :       ALLOCATE (helium%rho_property(rho_projected_area)%filename_suffix(nc))
    2376          30 :       ALLOCATE (helium%rho_property(rho_projected_area)%component_name(nc))
    2377          30 :       ALLOCATE (helium%rho_property(rho_projected_area)%component_index(nc))
    2378          30 :       helium%rho_property(rho_projected_area)%filename_suffix(1) = 'pa_x'
    2379          30 :       helium%rho_property(rho_projected_area)%filename_suffix(2) = 'pa_y'
    2380          30 :       helium%rho_property(rho_projected_area)%filename_suffix(3) = 'pa_z'
    2381          30 :       helium%rho_property(rho_projected_area)%component_name(1) = 'component x'
    2382          30 :       helium%rho_property(rho_projected_area)%component_name(2) = 'component y'
    2383          30 :       helium%rho_property(rho_projected_area)%component_name(3) = 'component z'
    2384         120 :       helium%rho_property(rho_projected_area)%component_index(:) = 0
    2385             : 
    2386          30 :       helium%rho_property(rho_winding_number)%name = 'Winding number squared density, W*W(r)'
    2387          30 :       nc = 3
    2388          30 :       helium%rho_property(rho_winding_number)%num_components = nc
    2389          30 :       ALLOCATE (helium%rho_property(rho_winding_number)%filename_suffix(nc))
    2390          30 :       ALLOCATE (helium%rho_property(rho_winding_number)%component_name(nc))
    2391          30 :       ALLOCATE (helium%rho_property(rho_winding_number)%component_index(nc))
    2392          30 :       helium%rho_property(rho_winding_number)%filename_suffix(1) = 'wn_x'
    2393          30 :       helium%rho_property(rho_winding_number)%filename_suffix(2) = 'wn_y'
    2394          30 :       helium%rho_property(rho_winding_number)%filename_suffix(3) = 'wn_z'
    2395          30 :       helium%rho_property(rho_winding_number)%component_name(1) = 'component x'
    2396          30 :       helium%rho_property(rho_winding_number)%component_name(2) = 'component y'
    2397          30 :       helium%rho_property(rho_winding_number)%component_name(3) = 'component z'
    2398         120 :       helium%rho_property(rho_winding_number)%component_index(:) = 0
    2399             : 
    2400          30 :       helium%rho_property(rho_winding_cycle)%name = 'Winding number squared density, W^2(r)'
    2401          30 :       nc = 3
    2402          30 :       helium%rho_property(rho_winding_cycle)%num_components = nc
    2403          30 :       ALLOCATE (helium%rho_property(rho_winding_cycle)%filename_suffix(nc))
    2404          30 :       ALLOCATE (helium%rho_property(rho_winding_cycle)%component_name(nc))
    2405          30 :       ALLOCATE (helium%rho_property(rho_winding_cycle)%component_index(nc))
    2406          30 :       helium%rho_property(rho_winding_cycle)%filename_suffix(1) = 'wc_x'
    2407          30 :       helium%rho_property(rho_winding_cycle)%filename_suffix(2) = 'wc_y'
    2408          30 :       helium%rho_property(rho_winding_cycle)%filename_suffix(3) = 'wc_z'
    2409          30 :       helium%rho_property(rho_winding_cycle)%component_name(1) = 'component x'
    2410          30 :       helium%rho_property(rho_winding_cycle)%component_name(2) = 'component y'
    2411          30 :       helium%rho_property(rho_winding_cycle)%component_name(3) = 'component z'
    2412         120 :       helium%rho_property(rho_winding_cycle)%component_index(:) = 0
    2413             : 
    2414          30 :       helium%rho_property(rho_moment_of_inertia)%name = 'Moment of inertia'
    2415          30 :       nc = 3
    2416          30 :       helium%rho_property(rho_moment_of_inertia)%num_components = nc
    2417          30 :       ALLOCATE (helium%rho_property(rho_moment_of_inertia)%filename_suffix(nc))
    2418          30 :       ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_name(nc))
    2419          30 :       ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_index(nc))
    2420          30 :       helium%rho_property(rho_moment_of_inertia)%filename_suffix(1) = 'mi_x'
    2421          30 :       helium%rho_property(rho_moment_of_inertia)%filename_suffix(2) = 'mi_y'
    2422          30 :       helium%rho_property(rho_moment_of_inertia)%filename_suffix(3) = 'mi_z'
    2423          30 :       helium%rho_property(rho_moment_of_inertia)%component_name(1) = 'component x'
    2424          30 :       helium%rho_property(rho_moment_of_inertia)%component_name(2) = 'component y'
    2425          30 :       helium%rho_property(rho_moment_of_inertia)%component_name(3) = 'component z'
    2426         120 :       helium%rho_property(rho_moment_of_inertia)%component_index(:) = 0
    2427             : 
    2428         180 :       helium%rho_property(:)%is_calculated = .FALSE.
    2429             : 
    2430          30 :       RETURN
    2431             :    END SUBROUTINE helium_rho_property_init
    2432             : 
    2433             : ! ***************************************************************************
    2434             : !> \brief  Create the density related data structures and initialize
    2435             : !> \param helium ...
    2436             : !> \date   2014-02-25
    2437             : !> \author Lukasz Walewski
    2438             : ! **************************************************************************************************
    2439          30 :    SUBROUTINE helium_rho_init(helium)
    2440             : 
    2441             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2442             : 
    2443             :       INTEGER                                            :: ii, itmp, jtmp
    2444             :       LOGICAL                                            :: explicit, ltmp
    2445             : 
    2446          30 :       CALL helium_rho_property_init(helium)
    2447             : 
    2448          30 :       helium%rho_num_act = 0
    2449             : 
    2450             :       ! check for atom number density
    2451             :       CALL section_vals_val_get( &
    2452             :          helium%input, &
    2453             :          "MOTION%PINT%HELIUM%RHO%ATOM_NUMBER", &
    2454          30 :          l_val=ltmp)
    2455          30 :       IF (ltmp) THEN
    2456          30 :          helium%rho_property(rho_atom_number)%is_calculated = .TRUE.
    2457          30 :          helium%rho_num_act = helium%rho_num_act + 1
    2458          30 :          helium%rho_property(rho_atom_number)%component_index(1) = helium%rho_num_act
    2459             :       END IF
    2460             : 
    2461             :       ! check for projected area density
    2462             :       CALL section_vals_val_get( &
    2463             :          helium%input, &
    2464             :          "MOTION%PINT%HELIUM%RHO%PROJECTED_AREA_2", &
    2465          30 :          l_val=ltmp)
    2466          30 :       IF (ltmp) THEN
    2467           0 :          helium%rho_property(rho_projected_area)%is_calculated = .TRUE.
    2468           0 :          DO ii = 1, helium%rho_property(rho_projected_area)%num_components
    2469           0 :             helium%rho_num_act = helium%rho_num_act + 1
    2470           0 :             helium%rho_property(rho_projected_area)%component_index(ii) = helium%rho_num_act
    2471             :          END DO
    2472             :       END IF
    2473             : 
    2474             :       ! check for winding number density
    2475             :       CALL section_vals_val_get( &
    2476             :          helium%input, &
    2477             :          "MOTION%PINT%HELIUM%RHO%WINDING_NUMBER_2", &
    2478          30 :          l_val=ltmp)
    2479          30 :       IF (ltmp) THEN
    2480           0 :          helium%rho_property(rho_winding_number)%is_calculated = .TRUE.
    2481           0 :          DO ii = 1, helium%rho_property(rho_winding_number)%num_components
    2482           0 :             helium%rho_num_act = helium%rho_num_act + 1
    2483           0 :             helium%rho_property(rho_winding_number)%component_index(ii) = helium%rho_num_act
    2484             :          END DO
    2485             :       END IF
    2486             : 
    2487             :       ! check for winding cycle density
    2488             :       CALL section_vals_val_get( &
    2489             :          helium%input, &
    2490             :          "MOTION%PINT%HELIUM%RHO%WINDING_CYCLE_2", &
    2491          30 :          l_val=ltmp)
    2492          30 :       IF (ltmp) THEN
    2493           0 :          helium%rho_property(rho_winding_cycle)%is_calculated = .TRUE.
    2494           0 :          DO ii = 1, helium%rho_property(rho_winding_cycle)%num_components
    2495           0 :             helium%rho_num_act = helium%rho_num_act + 1
    2496           0 :             helium%rho_property(rho_winding_cycle)%component_index(ii) = helium%rho_num_act
    2497             :          END DO
    2498             :       END IF
    2499             : 
    2500             :       ! check for moment of inertia density
    2501             :       CALL section_vals_val_get( &
    2502             :          helium%input, &
    2503             :          "MOTION%PINT%HELIUM%RHO%MOMENT_OF_INERTIA", &
    2504          30 :          l_val=ltmp)
    2505          30 :       IF (ltmp) THEN
    2506           0 :          helium%rho_property(rho_moment_of_inertia)%is_calculated = .TRUE.
    2507           0 :          DO ii = 1, helium%rho_property(rho_moment_of_inertia)%num_components
    2508           0 :             helium%rho_num_act = helium%rho_num_act + 1
    2509           0 :             helium%rho_property(rho_moment_of_inertia)%component_index(ii) = helium%rho_num_act
    2510             :          END DO
    2511             :       END IF
    2512             : 
    2513             :       ! set the cube dimensions, etc (common to all estimators)
    2514          30 :       helium%rho_maxr = helium%cell_size
    2515             :       CALL section_vals_val_get( &
    2516             :          helium%input, &
    2517             :          "MOTION%PINT%HELIUM%RHO%NBIN", &
    2518          30 :          i_val=helium%rho_nbin)
    2519          30 :       helium%rho_delr = helium%rho_maxr/REAL(helium%rho_nbin, dp)
    2520             : 
    2521             :       ! check for optional estimators based on winding paths
    2522          30 :       helium%rho_num_min_len_wdg = 0
    2523             :       CALL section_vals_val_get( &
    2524             :          helium%input, &
    2525             :          "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_WDG", &
    2526          30 :          explicit=explicit)
    2527          30 :       IF (explicit) THEN
    2528           0 :          NULLIFY (helium%rho_min_len_wdg_vals)
    2529             :          CALL section_vals_val_get( &
    2530             :             helium%input, &
    2531             :             "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_WDG", &
    2532           0 :             i_vals=helium%rho_min_len_wdg_vals)
    2533           0 :          itmp = SIZE(helium%rho_min_len_wdg_vals)
    2534           0 :          IF (itmp .GT. 0) THEN
    2535           0 :             helium%rho_num_min_len_wdg = itmp
    2536           0 :             helium%rho_num_act = helium%rho_num_act + itmp
    2537             :          END IF
    2538             :       END IF
    2539             : 
    2540             :       ! check for optional estimators based on non-winding paths
    2541          30 :       helium%rho_num_min_len_non = 0
    2542             :       CALL section_vals_val_get( &
    2543             :          helium%input, &
    2544             :          "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_NON", &
    2545          30 :          explicit=explicit)
    2546          30 :       IF (explicit) THEN
    2547           0 :          NULLIFY (helium%rho_min_len_non_vals)
    2548             :          CALL section_vals_val_get( &
    2549             :             helium%input, &
    2550             :             "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_NON", &
    2551           0 :             i_vals=helium%rho_min_len_non_vals)
    2552           0 :          itmp = SIZE(helium%rho_min_len_non_vals)
    2553           0 :          IF (itmp .GT. 0) THEN
    2554           0 :             helium%rho_num_min_len_non = itmp
    2555           0 :             helium%rho_num_act = helium%rho_num_act + itmp
    2556             :          END IF
    2557             :       END IF
    2558             : 
    2559             :       ! check for optional estimators based on all paths
    2560          30 :       helium%rho_num_min_len_all = 0
    2561             :       CALL section_vals_val_get( &
    2562             :          helium%input, &
    2563             :          "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_ALL", &
    2564          30 :          explicit=explicit)
    2565          30 :       IF (explicit) THEN
    2566           0 :          NULLIFY (helium%rho_min_len_all_vals)
    2567             :          CALL section_vals_val_get( &
    2568             :             helium%input, &
    2569             :             "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_ALL", &
    2570           0 :             i_vals=helium%rho_min_len_all_vals)
    2571           0 :          itmp = SIZE(helium%rho_min_len_all_vals)
    2572           0 :          IF (itmp .GT. 0) THEN
    2573           0 :             helium%rho_num_min_len_all = itmp
    2574           0 :             helium%rho_num_act = helium%rho_num_act + itmp
    2575             :          END IF
    2576             :       END IF
    2577             : 
    2578             :       ! get the weighting factor
    2579             :       CALL section_vals_val_get( &
    2580             :          helium%input, &
    2581             :          "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
    2582          30 :          i_val=helium%rho_iweight)
    2583             : 
    2584             :       ! allocate and initialize memory for density storage
    2585          30 :       itmp = helium%rho_nbin
    2586          30 :       jtmp = helium%rho_num_act
    2587         180 :       ALLOCATE (helium%rho_inst(jtmp, itmp, itmp, itmp))
    2588         180 :       ALLOCATE (helium%rho_accu(jtmp, itmp, itmp, itmp))
    2589         180 :       ALLOCATE (helium%rho_rstr(jtmp, itmp, itmp, itmp))
    2590         150 :       ALLOCATE (helium%rho_incr(jtmp, helium%atoms, helium%beads))
    2591             : 
    2592       28440 :       helium%rho_incr(:, :, :) = 0.0_dp
    2593    56287050 :       helium%rho_inst(:, :, :, :) = 0.0_dp
    2594    56287050 :       helium%rho_accu(:, :, :, :) = 0.0_dp
    2595    56287050 :       helium%rho_rstr(:, :, :, :) = 0.0_dp
    2596             : 
    2597          30 :       RETURN
    2598         240 :    END SUBROUTINE helium_rho_init
    2599             : 
    2600             : ! ***************************************************************************
    2601             : !> \brief  Restore the densities from the input structure.
    2602             : !> \param helium_env ...
    2603             : !> \date   2011-06-22
    2604             : !> \par    History
    2605             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    2606             : !> \author Lukasz Walewski
    2607             : ! **************************************************************************************************
    2608          24 :    SUBROUTINE helium_rho_restore(helium_env)
    2609             : 
    2610             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    2611             : 
    2612             :       CHARACTER(len=default_string_length)               :: stmp1, stmp2
    2613             :       CHARACTER(len=max_line_length)                     :: err_str
    2614             :       INTEGER                                            :: itmp, k, msglen
    2615             :       LOGICAL                                            :: explicit, ltmp
    2616          24 :       LOGICAL, DIMENSION(:, :, :, :), POINTER            :: m
    2617          24 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: message
    2618          24 :       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: f
    2619             : 
    2620             :       CALL section_vals_val_get(helium_env(1)%helium%input, &
    2621             :                                 "MOTION%PINT%HELIUM%AVERAGES%RHO", &
    2622          24 :                                 explicit=explicit)
    2623          24 :       IF (explicit) THEN
    2624           0 :          NULLIFY (message)
    2625             :          CALL section_vals_val_get(helium_env(1)%helium%input, &
    2626             :                                    "MOTION%PINT%HELIUM%AVERAGES%RHO", &
    2627           0 :                                    r_vals=message)
    2628           0 :          msglen = SIZE(message)
    2629           0 :          itmp = SIZE(helium_env(1)%helium%rho_rstr)
    2630           0 :          ltmp = (msglen .EQ. itmp)
    2631           0 :          IF (.NOT. ltmp) THEN
    2632           0 :             stmp1 = ""
    2633           0 :             WRITE (stmp1, *) msglen
    2634           0 :             stmp2 = ""
    2635           0 :             WRITE (stmp2, *) itmp
    2636             :             err_str = "Size of the S density array in the input ("// &
    2637             :                       TRIM(ADJUSTL(stmp1))// &
    2638             :                       ") .NE. that in the runtime ("// &
    2639           0 :                       TRIM(ADJUSTL(stmp2))//")."
    2640           0 :             CPABORT(err_str)
    2641             :          END IF
    2642             :       ELSE
    2643             :          RETURN
    2644             :       END IF
    2645             : 
    2646           0 :       itmp = helium_env(1)%helium%rho_nbin
    2647             :       NULLIFY (m, f)
    2648           0 :       ALLOCATE (m(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp))
    2649           0 :       ALLOCATE (f(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp))
    2650           0 :       m(:, :, :, :) = .TRUE.
    2651           0 :       f(:, :, :, :) = 0.0_dp
    2652             : 
    2653           0 :       DO k = 1, SIZE(helium_env)
    2654           0 :          helium_env(k)%helium%rho_rstr(:, :, :, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f)
    2655             :       END DO
    2656             : 
    2657           0 :       DEALLOCATE (f, m)
    2658             :       NULLIFY (message)
    2659             : 
    2660           0 :       RETURN
    2661          24 :    END SUBROUTINE helium_rho_restore
    2662             : 
    2663             : ! ***************************************************************************
    2664             : !> \brief Count atoms of different types and store their global indices.
    2665             : !> \param helium ...
    2666             : !> \param pint_env ...
    2667             : !> \author Lukasz Walewski
    2668             : !> \note  Arrays ALLOCATEd here are (should be) DEALLOCATEd in
    2669             : !>        helium_release.
    2670             : ! **************************************************************************************************
    2671          20 :    SUBROUTINE helium_set_solute_indices(helium, pint_env)
    2672             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2673             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2674             : 
    2675             :       INTEGER                                            :: iatom, natoms
    2676             :       REAL(KIND=dp)                                      :: mass
    2677             :       TYPE(cp_subsys_type), POINTER                      :: my_subsys
    2678             :       TYPE(f_env_type), POINTER                          :: my_f_env
    2679             :       TYPE(particle_list_type), POINTER                  :: my_particles
    2680             : 
    2681             : ! set up my_particles structure
    2682             : 
    2683          20 :       NULLIFY (my_f_env, my_subsys, my_particles)
    2684             :       CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
    2685          20 :                               f_env=my_f_env)
    2686          20 :       CALL force_env_get(force_env=my_f_env%force_env, subsys=my_subsys)
    2687          20 :       CALL cp_subsys_get(my_subsys, particles=my_particles)
    2688          20 :       CALL f_env_rm_defaults(my_f_env)
    2689             : 
    2690          20 :       natoms = helium%solute_atoms
    2691          20 :       NULLIFY (helium%solute_element)
    2692          40 :       ALLOCATE (helium%solute_element(natoms))
    2693             : 
    2694             :       ! find out how many different atomic types are there
    2695          20 :       helium%enum = 0
    2696          80 :       DO iatom = 1, natoms
    2697             :          CALL get_atomic_kind(my_particles%els(iatom)%atomic_kind, &
    2698             :                               mass=mass, &
    2699          80 :                               element_symbol=helium%solute_element(iatom))
    2700             :       END DO
    2701             : 
    2702          20 :       RETURN
    2703             :    END SUBROUTINE helium_set_solute_indices
    2704             : 
    2705             : ! ***************************************************************************
    2706             : !> \brief Sets helium%solute_cell based on the solute's force_env.
    2707             : !> \param helium ...
    2708             : !> \param pint_env ...
    2709             : !> \author Lukasz Walewski
    2710             : !> \note  The simulation cell for the solvated molecule is taken from force_env
    2711             : !>        which should assure that we get proper cell dimensions regardless of
    2712             : !>        the method used for the solute (QS, FIST). Helium solvent needs the
    2713             : !>        solute's cell dimensions to calculate the solute-solvent distances
    2714             : !>        correctly.
    2715             : !> \note  At the moment only orthorhombic cells are supported.
    2716             : ! **************************************************************************************************
    2717          40 :    SUBROUTINE helium_set_solute_cell(helium, pint_env)
    2718             :       TYPE(helium_solvent_type), POINTER                 :: helium
    2719             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
    2720             : 
    2721             :       LOGICAL                                            :: my_orthorhombic
    2722             :       TYPE(cell_type), POINTER                           :: my_cell
    2723             :       TYPE(f_env_type), POINTER                          :: my_f_env
    2724             : 
    2725             : ! get the cell structure from pint_env
    2726             : 
    2727          20 :       NULLIFY (my_f_env, my_cell)
    2728             :       CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
    2729          20 :                               f_env=my_f_env)
    2730          20 :       CALL force_env_get(force_env=my_f_env%force_env, cell=my_cell)
    2731          20 :       CALL f_env_rm_defaults(my_f_env)
    2732             : 
    2733          20 :       CALL get_cell(my_cell, orthorhombic=my_orthorhombic)
    2734          20 :       IF (.NOT. my_orthorhombic) THEN
    2735           0 :          CPABORT("Helium solvent not implemented for non-orthorhombic cells.")
    2736             :       ELSE
    2737          20 :          helium%solute_cell => my_cell
    2738             :       END IF
    2739             : 
    2740          20 :       RETURN
    2741             :    END SUBROUTINE helium_set_solute_cell
    2742             : 
    2743             : END MODULE helium_methods

Generated by: LCOV version 1.15