LCOV - code coverage report
Current view: top level - src/motion - helium_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 1097 1354 81.0 %
Date: 2024-12-12 07:20:04 Functions: 22 22 100.0 %

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

Generated by: LCOV version 1.15