LCOV - code coverage report
Current view: top level - src/motion - helium_interactions.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 273 332 82.2 %
Date: 2024-04-20 06:29:22 Functions: 10 11 90.9 %

          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 that handle helium-solvent and helium-helium interactions
      10             : !> \author Lukasz Walewski
      11             : !> \date   2009-06-10
      12             : ! **************************************************************************************************
      13             : MODULE helium_interactions
      14             : 
      15             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      16             :                                               cp_logger_type
      17             :    USE helium_common,                   ONLY: helium_eval_chain,&
      18             :                                               helium_eval_expansion,&
      19             :                                               helium_pbc,&
      20             :                                               helium_spline
      21             :    USE helium_nnp,                      ONLY: helium_nnp_print
      22             :    USE helium_types,                    ONLY: e_id_interact,&
      23             :                                               e_id_kinetic,&
      24             :                                               e_id_potential,&
      25             :                                               e_id_thermo,&
      26             :                                               e_id_total,&
      27             :                                               e_id_virial,&
      28             :                                               helium_solvent_p_type,&
      29             :                                               helium_solvent_type
      30             :    USE input_constants,                 ONLY: helium_sampling_worm,&
      31             :                                               helium_solute_intpot_mwater,&
      32             :                                               helium_solute_intpot_nnp,&
      33             :                                               helium_solute_intpot_none
      34             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      35             :                                               section_vals_type
      36             :    USE kinds,                           ONLY: dp
      37             :    USE nnp_acsf,                        ONLY: nnp_calc_acsf
      38             :    USE nnp_environment_types,           ONLY: nnp_type
      39             :    USE nnp_model,                       ONLY: nnp_gradients,&
      40             :                                               nnp_predict
      41             :    USE physcon,                         ONLY: angstrom,&
      42             :                                               kelvin
      43             :    USE pint_types,                      ONLY: pint_env_type
      44             :    USE splines_types,                   ONLY: spline_data_type
      45             : #include "../base/base_uses.f90"
      46             : 
      47             :    IMPLICIT NONE
      48             : 
      49             :    PRIVATE
      50             : 
      51             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      52             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_interactions'
      53             : 
      54             :    PUBLIC :: helium_calc_energy
      55             :    PUBLIC :: helium_total_link_action
      56             :    PUBLIC :: helium_total_pair_action
      57             :    PUBLIC :: helium_total_inter_action
      58             :    PUBLIC :: helium_solute_e_f
      59             :    PUBLIC :: helium_bead_solute_e_f
      60             :    PUBLIC :: helium_intpot_scan
      61             :    PUBLIC :: helium_vij
      62             : 
      63             : CONTAINS
      64             : 
      65             : ! ***************************************************************************
      66             : !> \brief  Calculate the helium energy (including helium-solute interaction)
      67             : !> \param    helium     helium environment
      68             : !> \param    pint_env   path integral environment
      69             : !> \par History
      70             : !>         2009-06 moved I/O out from here [lwalewski]
      71             : !> \author hforbert
      72             : ! **************************************************************************************************
      73        7047 :    SUBROUTINE helium_calc_energy(helium, pint_env)
      74             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
      75             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
      76             : 
      77             :       INTEGER                                            :: b, bead, i, j, n
      78        7047 :       INTEGER, DIMENSION(:), POINTER                     :: perm
      79             :       LOGICAL                                            :: nperiodic
      80             :       REAL(KIND=dp)                                      :: a, cell_size, en, interac, kin, pot, &
      81             :                                                             rmax, rmin, vkin
      82        7047 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work2, work3
      83        7047 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
      84             :       REAL(KIND=dp), DIMENSION(3)                        :: r
      85        7047 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos
      86             :       TYPE(spline_data_type), POINTER                    :: e0
      87             : 
      88        7047 :       pos => helium%pos
      89        7047 :       perm => helium%permutation
      90        7047 :       e0 => helium%e0
      91        7047 :       cell_size = 0.5_dp*helium%cell_size
      92        7047 :       nperiodic = .NOT. helium%periodic
      93        7047 :       n = helium%atoms
      94        7047 :       b = helium%beads
      95        7047 :       en = 0.0_dp
      96        7047 :       pot = 0.0_dp
      97        7047 :       rmin = 1.0e20_dp
      98        7047 :       rmax = 0.0_dp
      99             :       ALLOCATE (work(3, helium%beads + 1), &
     100             :                 work2(helium%beads + 1), &
     101       49329 :                 work3(SIZE(helium%uoffdiag, 1) + 1))
     102      214104 :       DO i = 1, n - 1
     103     3503766 :          DO j = i + 1, n
     104    71146494 :             DO bead = 1, b
     105   274716990 :                work(:, bead) = pos(:, i, bead) - pos(:, j, bead)
     106             :             END DO
     107    13158648 :             work(:, b + 1) = pos(:, perm(i), 1) - pos(:, perm(j), 1)
     108     3289662 :             en = en + helium_eval_chain(helium, work, b + 1, work2, work3, energy=.TRUE.)
     109    71353551 :             DO bead = 1, b
     110    67856832 :                a = work2(bead)
     111    67856832 :                IF (a < rmin) rmin = a
     112    67856832 :                IF (a > rmax) rmax = a
     113    71146494 :                IF ((a < cell_size) .OR. nperiodic) THEN
     114    63532837 :                   pot = pot + helium_spline(helium%vij, a)
     115             :                END IF
     116             :             END DO
     117             :          END DO
     118             :       END DO
     119        7047 :       DEALLOCATE (work, work2, work3)
     120        7047 :       pot = pot/b
     121        7047 :       en = en/b
     122             : 
     123             :       ! helium-solute interaction energy (all beads of all particles)
     124        7047 :       interac = 0.0_dp
     125        7047 :       IF (helium%solute_present) THEN
     126        3637 :          CALL helium_solute_e(pint_env, helium, interac)
     127             :       END IF
     128        7047 :       interac = interac/b
     129             : 
     130             : !TODO:
     131        7047 :       vkin = 0.0_dp
     132             : !   vkin = helium_virial_energy(helium)
     133             : 
     134        7047 :       kin = 0.0_dp
     135      221151 :       DO i = 1, n
     136      856416 :          r(:) = pos(:, i, b) - pos(:, perm(i), 1)
     137      214104 :          CALL helium_pbc(helium, r)
     138      214104 :          kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     139     4414791 :          DO bead = 2, b
     140    16774560 :             r(:) = pos(:, i, bead - 1) - pos(:, i, bead)
     141     4193640 :             CALL helium_pbc(helium, r)
     142     4407744 :             kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     143             :          END DO
     144             :       END DO
     145        7047 :       kin = 1.5_dp*n/helium%tau - 0.5*kin/(b*helium%tau**2*helium%hb2m)
     146             : 
     147             : ! TODO: move printing somewhere else ?
     148             : !   print *,"POT = ",(pot/n+helium%e_corr)*kelvin,"K"
     149             : !   print *,"INTERAC = ",interac*kelvin,"K"
     150             : !   print *,"RMIN= ",rmin*angstrom,"A"
     151             : !   print *,"RMAX= ",rmax*angstrom,"A"
     152             : !   print *,"EVIRIAL not valid!"
     153             : !   print *,"ETHERMO= ",((en+kin)/n+helium%e_corr)*kelvin,"K"
     154             : !   print *,"ECORR= ",helium%e_corr*kelvin,"K"
     155             : !!   kin = helium_total_action(helium)
     156             : !!   print *,"ACTION= ",kin
     157             : !   print *,"WINDING#= ",helium_calc_winding(helium)
     158             : 
     159        7047 :       helium%energy_inst(e_id_potential) = pot/n + helium%e_corr
     160        7047 :       helium%energy_inst(e_id_kinetic) = (en - pot + kin)/n
     161        7047 :       helium%energy_inst(e_id_interact) = interac
     162        7047 :       helium%energy_inst(e_id_thermo) = (en + kin)/n + helium%e_corr
     163        7047 :       helium%energy_inst(e_id_virial) = vkin ! 0.0_dp at the moment
     164        7047 :       helium%energy_inst(e_id_total) = helium%energy_inst(e_id_thermo)
     165             :       ! Once vkin is properly implemented, switch to:
     166             :       ! helium%energy_inst(e_id_total) = (en+vkin)/n+helium%e_corr
     167             : 
     168       14094 :    END SUBROUTINE helium_calc_energy
     169             : 
     170             : ! ***************************************************************************
     171             : !> \brief  Computes the total harmonic link action of the helium
     172             : !> \param helium ...
     173             : !> \return ...
     174             : !> \date   2016-05-03
     175             : !> \author Felix Uhl
     176             : ! **************************************************************************************************
     177          50 :    REAL(KIND=dp) FUNCTION helium_total_link_action(helium) RESULT(linkaction)
     178             : 
     179             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     180             : 
     181             :       INTEGER                                            :: iatom, ibead
     182          50 :       INTEGER, DIMENSION(:), POINTER                     :: perm
     183             :       REAL(KIND=dp), DIMENSION(3)                        :: r
     184             : 
     185          50 :       perm => helium%permutation
     186          50 :       linkaction = 0.0_dp
     187             : 
     188             :       ! Harmonic Link action
     189             :       ! (r(m-1) - r(m))**2/(4*lambda*tau)
     190         800 :       DO ibead = 1, helium%beads - 1
     191        4550 :          DO iatom = 1, helium%atoms
     192       15000 :             r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, iatom, ibead + 1)
     193        3750 :             CALL helium_pbc(helium, r)
     194        4500 :             linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     195             :          END DO
     196             :       END DO
     197         300 :       DO iatom = 1, helium%atoms
     198             :          ! choose last bead connection according to permutation table
     199        1000 :          r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, perm(iatom), 1)
     200         250 :          CALL helium_pbc(helium, r)
     201         300 :          linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     202             :       END DO
     203          50 :       linkaction = linkaction/(2.0_dp*helium%tau*helium%hb2m)
     204             : 
     205          50 :    END FUNCTION helium_total_link_action
     206             : 
     207             : ! ***************************************************************************
     208             : !> \brief  Computes the total pair action of the helium
     209             : !> \param helium ...
     210             : !> \return ...
     211             : !> \date   2016-05-03
     212             : !> \author Felix Uhl
     213             : ! **************************************************************************************************
     214          50 :    REAL(KIND=dp) FUNCTION helium_total_pair_action(helium) RESULT(pairaction)
     215             : 
     216             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     217             : 
     218             :       INTEGER                                            :: iatom, ibead, jatom, opatom, patom
     219          50 :       INTEGER, DIMENSION(:), POINTER                     :: perm
     220          50 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     221             :       REAL(KIND=dp), DIMENSION(3)                        :: r, rp
     222             : 
     223         150 :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
     224          50 :       perm => helium%permutation
     225          50 :       pairaction = 0.0_dp
     226             : 
     227             :       ! He-He pair action
     228         800 :       DO ibead = 1, helium%beads - 1
     229        3800 :          DO iatom = 1, helium%atoms - 1
     230       11250 :             DO jatom = iatom + 1, helium%atoms
     231       30000 :                r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
     232       30000 :                rp(:) = helium%pos(:, iatom, ibead + 1) - helium%pos(:, jatom, ibead + 1)
     233       10500 :                pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     234             :             END DO
     235             :          END DO
     236             :       END DO
     237             :       !Ensure right permutation for pair action of last and first beads.
     238         250 :       DO iatom = 1, helium%atoms - 1
     239         750 :          DO jatom = iatom + 1, helium%atoms
     240        2000 :             r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, jatom, helium%beads)
     241        2000 :             rp(:) = helium%pos(:, perm(iatom), 1) - helium%pos(:, perm(jatom), 1)
     242         700 :             pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     243             :          END DO
     244             :       END DO
     245             : 
     246             :       ! correct for open worm configurations
     247          50 :       IF (.NOT. helium%worm_is_closed) THEN
     248             :          ! special treatment if double bead is first bead
     249           0 :          iatom = helium%worm_atom_idx
     250           0 :          IF (helium%worm_bead_idx == 1) THEN
     251             :             ! patom is the atom in front of the lone head bead
     252           0 :             patom = helium%iperm(iatom)
     253             :             ! go through all atoms
     254           0 :             DO jatom = 1, helium%atoms
     255           0 :                IF (jatom == helium%worm_atom_idx) CYCLE
     256           0 :                opatom = helium%iperm(jatom)
     257             :                ! subtract pair action for closed link
     258           0 :                r(:) = helium%pos(:, iatom, 1) - helium%pos(:, jatom, 1)
     259           0 :                rp(:) = helium%pos(:, patom, helium%beads) - helium%pos(:, opatom, helium%beads)
     260           0 :                pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
     261             :                ! and add corrected extra link
     262             :                ! rp stays the same
     263           0 :                r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, 1)
     264           0 :                pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     265             :             END DO
     266             :          ELSE
     267             :             ! bead stays constant
     268           0 :             ibead = helium%worm_bead_idx
     269             :             ! go through all atoms
     270           0 :             DO jatom = 1, helium%atoms
     271           0 :                IF (jatom == helium%worm_atom_idx) CYCLE
     272             :                ! subtract pair action for closed link
     273           0 :                r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
     274           0 :                rp(:) = helium%pos(:, iatom, ibead - 1) - helium%pos(:, jatom, ibead - 1)
     275           0 :                pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
     276             :                ! and add corrected extra link
     277             :                ! rp stays the same
     278           0 :                r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, ibead)
     279           0 :                pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     280             :             END DO
     281             :          END IF
     282             :       END IF
     283          50 :       DEALLOCATE (work3)
     284             : 
     285          50 :    END FUNCTION helium_total_pair_action
     286             : 
     287             : ! ***************************************************************************
     288             : !> \brief  Computes the total interaction of the helium with the solute
     289             : !> \param pint_env ...
     290             : !> \param helium ...
     291             : !> \return ...
     292             : !> \date   2016-05-03
     293             : !> \author Felix Uhl
     294             : ! **************************************************************************************************
     295          50 :    REAL(KIND=dp) FUNCTION helium_total_inter_action(pint_env, helium) RESULT(interaction)
     296             : 
     297             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     298             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     299             : 
     300             :       INTEGER                                            :: iatom, ibead
     301             :       REAL(KIND=dp)                                      :: e
     302             : 
     303          50 :       interaction = 0.0_dp
     304             : 
     305             :       ! InterAction with solute
     306          50 :       IF (helium%solute_present) THEN
     307         850 :          DO ibead = 1, helium%beads
     308        4850 :             DO iatom = 1, helium%atoms
     309             : 
     310             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     311        4000 :                                            iatom, ibead, helium%pos(:, iatom, ibead), e)
     312        4800 :                interaction = interaction + e
     313             :             END DO
     314             :          END DO
     315          50 :          IF (helium%sampling_method == helium_sampling_worm) THEN
     316           0 :             IF (.NOT. helium%worm_is_closed) THEN
     317             :                ! subtract half of tail bead interaction again
     318             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     319             :                                            helium%worm_atom_idx, helium%worm_bead_idx, &
     320           0 :                                            helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), e)
     321           0 :                interaction = interaction - 0.5_dp*e
     322             :                ! add half of head bead interaction
     323             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     324             :                                            helium%worm_atom_idx, helium%worm_bead_idx, &
     325           0 :                                            helium%worm_xtra_bead, e)
     326           0 :                interaction = interaction + 0.5_dp*e
     327             :             END IF
     328             :          END IF
     329             :       END IF
     330             : 
     331          50 :       interaction = interaction*helium%tau
     332             : 
     333          50 :    END FUNCTION helium_total_inter_action
     334             : 
     335             : ! ***************************************************************************
     336             : !> \brief Calculate general helium-solute interaction energy (and forces)
     337             : !>        between one helium bead and the corresponding solute time slice.
     338             : !> \param pint_env           path integral environment
     339             : !> \param helium ...
     340             : !> \param helium_part_index  helium particle index
     341             : !> \param helium_slice_index helium time slice index
     342             : !> \param helium_r_opt       explicit helium bead coordinates (optional)
     343             : !> \param energy             calculated energy
     344             : !> \param force              calculated force (if requested)
     345             : !> \par History
     346             : !>         2019-09 Added multiple-time striding in imag. time [cschran]
     347             : !>         2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
     348             : !> \author Lukasz Walewski
     349             : ! **************************************************************************************************
     350     5542964 :    SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
     351             :                                      helium_slice_index, helium_r_opt, energy, force)
     352             : 
     353             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     354             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     355             :       INTEGER, INTENT(IN)                                :: helium_part_index, helium_slice_index
     356             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL  :: helium_r_opt
     357             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     358             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     359             :          OPTIONAL, POINTER                               :: force
     360             : 
     361             :       INTEGER                                            :: hbeads, hi, qi, stride
     362             :       REAL(KIND=dp), DIMENSION(3)                        :: helium_r
     363     5542964 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: my_force
     364             : 
     365     5542964 :       hbeads = helium%beads
     366             :       ! helium bead index that is invariant wrt the rotations
     367     5542964 :       hi = MOD(helium_slice_index - 1 + hbeads + helium%relrot, hbeads) + 1
     368             :       ! solute bead index that belongs to hi helium index
     369     5542964 :       qi = ((hi - 1)*pint_env%p)/hbeads + 1
     370             : 
     371             :       ! coordinates of the helium bead
     372     5542964 :       IF (PRESENT(helium_r_opt)) THEN
     373     1460066 :          helium_r(:) = helium_r_opt(:)
     374             :       ELSE
     375    16331592 :          helium_r(:) = helium%pos(:, helium_part_index, helium_slice_index)
     376             :       END IF
     377             : 
     378    11030756 :       SELECT CASE (helium%solute_interaction)
     379             : 
     380             :       CASE (helium_solute_intpot_mwater)
     381     5487792 :          IF (PRESENT(force)) THEN
     382    54652288 :             force(:, :) = 0.0_dp
     383     1171264 :             my_force => force(qi, :)
     384             :             CALL helium_intpot_model_water( &
     385             :                pint_env%x(qi, :), &
     386             :                helium, &
     387             :                helium_r, &
     388             :                energy, &
     389             :                my_force &
     390     1171264 :                )
     391             :          ELSE
     392             :             CALL helium_intpot_model_water( &
     393             :                pint_env%x(qi, :), &
     394             :                helium, &
     395             :                helium_r, &
     396             :                energy &
     397     4316528 :                )
     398             :          END IF
     399             : 
     400             :       CASE (helium_solute_intpot_nnp)
     401       55172 :          IF (PRESENT(force)) THEN
     402      246400 :             force(:, :) = 0.0_dp
     403        1600 :             my_force => force(qi, :)
     404             :             CALL helium_intpot_nnp( &
     405             :                pint_env%x(qi, :), &
     406             :                helium, &
     407             :                helium_r, &
     408             :                energy, &
     409             :                my_force &
     410        1600 :                )
     411             :          ELSE
     412             :             CALL helium_intpot_nnp( &
     413             :                pint_env%x(qi, :), &
     414             :                helium, &
     415             :                helium_r, &
     416             :                energy &
     417       53572 :                )
     418             :          END IF
     419             : 
     420             :       CASE (helium_solute_intpot_none)
     421           0 :          energy = 0.0_dp
     422     5542964 :          IF (PRESENT(force)) THEN
     423           0 :             force(:, :) = 0.0_dp
     424             :          END IF
     425             : 
     426             :       CASE DEFAULT
     427             : 
     428             :       END SELECT
     429             : 
     430             :       ! Account for Imaginary time striding in forces:
     431     5542964 :       IF (PRESENT(force)) THEN
     432     1172864 :          IF (hbeads < pint_env%p) THEN
     433        3072 :             stride = pint_env%p/hbeads
     434      915456 :             force = force*REAL(stride, dp)
     435             :          END IF
     436             :       END IF
     437             : 
     438     5542964 :    END SUBROUTINE helium_bead_solute_e_f
     439             : 
     440             : ! ***************************************************************************
     441             : !> \brief Calculate total helium-solute interaction energy and forces.
     442             : !> \param   pint_env   path integral environment
     443             : !> \param helium ...
     444             : !> \param   energy     calculated interaction energy
     445             : !> \author Lukasz Walewski
     446             : ! **************************************************************************************************
     447        2647 :    SUBROUTINE helium_solute_e_f(pint_env, helium, energy)
     448             : 
     449             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     450             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     451             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     452             : 
     453             :       INTEGER                                            :: ia, ib, jb, jc
     454             :       REAL(KIND=dp)                                      :: my_energy
     455        2647 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: force
     456             : 
     457        2647 :       NULLIFY (force)
     458        2647 :       force => helium%force_inst
     459             : 
     460        2647 :       energy = 0.0_dp
     461      123814 :       force(:, :) = 0.0_dp
     462             : 
     463             :       ! calculate the total interaction energy and gradients between the
     464             :       ! solute and the helium, sum over all beads of all He particles
     465       75951 :       DO ia = 1, helium%atoms
     466     1248815 :          DO ib = 1, helium%beads
     467             :             CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, &
     468     1172864 :                                         energy=my_energy, force=helium%rtmp_p_ndim_2d)
     469     1172864 :             energy = energy + my_energy
     470     6042840 :             DO jb = 1, pint_env%p
     471    49139584 :                DO jc = 1, pint_env%ndim
     472    47966720 :                   force(jb, jc) = force(jb, jc) + helium%rtmp_p_ndim_2d(jb, jc)
     473             :                END DO
     474             :             END DO
     475             :          END DO
     476             :       END DO
     477             : 
     478        2647 :    END SUBROUTINE helium_solute_e_f
     479             : 
     480             : ! ***************************************************************************
     481             : !> \brief Calculate total helium-solute interaction energy.
     482             : !> \param   pint_env   path integral environment
     483             : !> \param helium ...
     484             : !> \param   energy     calculated interaction energy
     485             : !> \author Lukasz Walewski
     486             : ! **************************************************************************************************
     487        3637 :    SUBROUTINE helium_solute_e(pint_env, helium, energy)
     488             : 
     489             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     490             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     491             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     492             : 
     493             :       INTEGER                                            :: ia, ib
     494             :       REAL(KIND=dp)                                      :: my_energy
     495             : 
     496        3637 :       energy = 0.0_dp
     497             : 
     498      108621 :       DO ia = 1, helium%atoms
     499     1788365 :          DO ib = 1, helium%beads
     500     1679744 :             CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, energy=my_energy)
     501     1784728 :             energy = energy + my_energy
     502             :          END DO
     503             :       END DO
     504             : 
     505        3637 :    END SUBROUTINE helium_solute_e
     506             : 
     507             : ! ***************************************************************************
     508             : !> \brief  Scan the helium-solute interaction energy within the periodic cell
     509             : !> \param pint_env ...
     510             : !> \param helium_env ...
     511             : !> \date   2014-01-22
     512             : !> \par    History
     513             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     514             : !> \author Lukasz Walewski
     515             : ! **************************************************************************************************
     516           0 :    SUBROUTINE helium_intpot_scan(pint_env, helium_env)
     517             : 
     518             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     519             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     520             : 
     521             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_intpot_scan'
     522             : 
     523             :       INTEGER                                            :: handle, ic, ix, iy, iz, k, nbin
     524             :       LOGICAL                                            :: wrapped
     525             :       REAL(KIND=dp)                                      :: delr, my_en, ox, oy, oz
     526             :       REAL(kind=dp), DIMENSION(3)                        :: pbc1, pbc2, pos
     527             : 
     528           0 :       CALL timeset(routineN, handle)
     529             : 
     530             :       ! Perform scan only on ionode, since this is only used to output the intpot
     531           0 :       IF (pint_env%logger%para_env%is_source()) THEN
     532             :          ! Assume ionode always to have at least one helium_env
     533           0 :          k = 1
     534           0 :          helium_env(k)%helium%rho_inst(1, :, :, :) = 0.0_dp
     535           0 :          nbin = helium_env(k)%helium%rho_nbin
     536           0 :          delr = helium_env(k)%helium%rho_delr
     537           0 :          helium_env(k)%helium%center(:) = 0.0_dp
     538           0 :          ox = helium_env(k)%helium%center(1) - helium_env(k)%helium%rho_maxr/2.0_dp
     539           0 :          oy = helium_env(k)%helium%center(2) - helium_env(k)%helium%rho_maxr/2.0_dp
     540           0 :          oz = helium_env(k)%helium%center(3) - helium_env(k)%helium%rho_maxr/2.0_dp
     541             : 
     542           0 :          DO ix = 1, nbin
     543           0 :             DO iy = 1, nbin
     544           0 :                DO iz = 1, nbin
     545             : 
     546             :                   ! put the probe in the center of the current voxel
     547           0 :                   pos(:) = (/ox + (ix - 0.5_dp)*delr, oy + (iy - 0.5_dp)*delr, oz + (iz - 0.5_dp)*delr/)
     548             : 
     549             :                   ! calc interaction energy for the current probe position
     550           0 :                   helium_env(k)%helium%pos(:, 1, 1) = pos(:)
     551           0 :                   CALL helium_bead_solute_e_f(pint_env, helium_env(k)%helium, 1, 1, energy=my_en)
     552             : 
     553             :                   ! check if the probe fits within the unit cell
     554           0 :                   pbc1(:) = pos(:) - helium_env(k)%helium%center
     555           0 :                   pbc2(:) = pbc1(:)
     556           0 :                   CALL helium_pbc(helium_env(k)%helium, pbc2)
     557           0 :                   wrapped = .FALSE.
     558           0 :                   DO ic = 1, 3
     559           0 :                      IF (ABS(pbc1(ic) - pbc2(ic)) .GT. 10.0_dp*EPSILON(0.0_dp)) THEN
     560           0 :                         wrapped = .TRUE.
     561             :                      END IF
     562             :                   END DO
     563             : 
     564             :                   ! set the interaction energy value
     565           0 :                   IF (wrapped) THEN
     566           0 :                      helium_env(k)%helium%rho_inst(1, ix, iy, iz) = 0.0_dp
     567             :                   ELSE
     568           0 :                      helium_env(k)%helium%rho_inst(1, ix, iy, iz) = my_en
     569             :                   END IF
     570             : 
     571             :                END DO
     572             :             END DO
     573             :          END DO
     574             :       END IF
     575             : 
     576           0 :       CALL timestop(handle)
     577           0 :    END SUBROUTINE helium_intpot_scan
     578             : 
     579             : ! ***************************************************************************
     580             : !> \brief Calculate model helium-solute interaction energy and forces
     581             : !>        between one helium bead and the corresponding solute time
     582             : !>        slice asuming water solute.
     583             : !> \param solute_x  solute positions ARR(3*NATOMS)
     584             : !>        to global atom indices
     585             : !> \param helium    only needed for helium_pbc call at the moment
     586             : !> \param helium_x  helium bead position ARR(3)
     587             : !> \param energy    calculated interaction energy
     588             : !> \param force ...
     589             : !> \author Felix Uhl
     590             : ! **************************************************************************************************
     591     5487792 :    SUBROUTINE helium_intpot_model_water(solute_x, helium, helium_x, energy, force)
     592             : 
     593             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: solute_x
     594             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     595             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: helium_x
     596             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     597             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
     598             :          OPTIONAL, POINTER                               :: force
     599             : 
     600             :       INTEGER                                            :: i, ig
     601             :       REAL(KIND=dp)                                      :: d, d2, dd, ep, eps, s1, s2, sig
     602             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, solute_r
     603             : 
     604     5487792 :       energy = 0.0_dp
     605     5487792 :       IF (PRESENT(force)) THEN
     606    11712640 :          force(:) = 0.0_dp
     607             :       END IF
     608             : 
     609     5487792 :       sig = 2.69_dp ! 1.4 Angstrom
     610     5487792 :       eps = 60.61e-6_dp ! 19 K
     611     5487792 :       s1 = 0.0_dp
     612    21951168 :       DO i = 1, SIZE(helium%solute_element)
     613    21951168 :          IF (helium%solute_element(i) == "H ") THEN
     614    10975584 :             ig = i - 1
     615    10975584 :             solute_r(1) = solute_x(3*ig + 1)
     616    10975584 :             solute_r(2) = solute_x(3*ig + 2)
     617    10975584 :             solute_r(3) = solute_x(3*ig + 3)
     618    43902336 :             dr(:) = solute_r(:) - helium_x(:)
     619    10975584 :             CALL helium_pbc(helium, dr)
     620    10975584 :             d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
     621    10975584 :             d = SQRT(d2)
     622    10975584 :             dd = (sig/d)**6
     623    10975584 :             ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
     624    10975584 :             s1 = s1 + ep
     625    10975584 :             s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
     626    10975584 :             IF (PRESENT(force)) THEN
     627     2342528 :                force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
     628     2342528 :                force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
     629     2342528 :                force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
     630             :             END IF
     631             :          END IF
     632             :       END DO ! i = 1, num_hydrogen
     633     5487792 :       energy = energy + s1
     634             : 
     635     5487792 :       sig = 5.01_dp ! 2.6 Angstrom
     636     5487792 :       eps = 104.5e-6_dp ! 33 K
     637     5487792 :       s1 = 0.0_dp
     638    21951168 :       DO i = 1, SIZE(helium%solute_element)
     639    21951168 :          IF (helium%solute_element(i) == "O ") THEN
     640     5487792 :             ig = i - 1
     641     5487792 :             solute_r(1) = solute_x(3*ig + 1)
     642     5487792 :             solute_r(2) = solute_x(3*ig + 2)
     643     5487792 :             solute_r(3) = solute_x(3*ig + 3)
     644    21951168 :             dr(:) = solute_r(:) - helium_x(:)
     645     5487792 :             CALL helium_pbc(helium, dr)
     646     5487792 :             d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
     647     5487792 :             d = SQRT(d2)
     648     5487792 :             dd = (sig/d)**6
     649     5487792 :             ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
     650     5487792 :             s1 = s1 + ep
     651     5487792 :             s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
     652     5487792 :             IF (PRESENT(force)) THEN
     653     1171264 :                force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
     654     1171264 :                force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
     655     1171264 :                force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
     656             :             END IF
     657             :          END IF
     658             :       END DO ! i = 1, num_chlorine
     659     5487792 :       energy = energy + s1
     660             : 
     661     5487792 :    END SUBROUTINE helium_intpot_model_water
     662             : 
     663             : ! ***************************************************************************
     664             : !> \brief  Calculate helium-solute interaction energy and forces between one
     665             : !>         helium bead and the corresponding solute time slice using NNP.
     666             : !> \param  solute_x  solute positions ARR(3*NATOMS)
     667             : !>         to global atom indices
     668             : !> \param  helium    only needed for helium_pbc call at the moment
     669             : !> \param  helium_x  helium bead position ARR(3)
     670             : !> \param  energy    calculated interaction energy
     671             : !> \param  force     (optional) calculated force
     672             : !> \date   2023-02-22
     673             : !> \author Laura Duran
     674             : ! **************************************************************************************************
     675       55172 :    SUBROUTINE helium_intpot_nnp(solute_x, helium, helium_x, energy, force)
     676             : 
     677             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: solute_x
     678             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     679             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: helium_x
     680             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     681             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
     682             :          OPTIONAL, POINTER                               :: force
     683             : 
     684             :       INTEGER                                            :: i, i_com, ig, ind, ind_he, j, k, m
     685             :       LOGICAL                                            :: extrapolate
     686             :       REAL(KIND=dp)                                      :: rsqr, rvect(3), threshold
     687       55172 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: denergydsym
     688       55172 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dsymdxyz
     689             :       TYPE(cp_logger_type), POINTER                      :: logger
     690             :       TYPE(nnp_type), POINTER                            :: nnp
     691             :       TYPE(section_vals_type), POINTER                   :: print_section
     692             : 
     693       55172 :       NULLIFY (logger)
     694       55172 :       logger => cp_get_default_logger()
     695             : 
     696       55172 :       IF (PRESENT(force)) THEN
     697       28800 :          helium%nnp%myforce(:, :, :) = 0.0_dp
     698             :       END IF
     699             : 
     700       55172 :       extrapolate = .FALSE.
     701       55172 :       threshold = 0.0001d0
     702             : 
     703             :       !fill coord array
     704       55172 :       ig = 1
     705      220688 :       DO i = 1, helium%nnp%n_ele
     706      165516 :          IF (helium%nnp%ele(i) == 'He') THEN
     707       55172 :             ind_he = ig
     708      220688 :             DO m = 1, 3
     709      220688 :                helium%nnp%coord(m, ig) = helium_x(m)
     710             :             END DO
     711       55172 :             ig = ig + 1
     712             :          END IF
     713      717236 :          DO j = 1, helium%solute_atoms
     714      662064 :             IF (helium%nnp%ele(i) == helium%solute_element(j)) THEN
     715      662064 :                DO m = 1, 3
     716      662064 :                   helium%nnp%coord(m, ig) = solute_x(3*(j - 1) + m)
     717             :                END DO
     718      165516 :                ig = ig + 1
     719             :             END IF
     720             :          END DO
     721             :       END DO
     722             : 
     723             :       ! check for hard core condition
     724       55172 :       IF (ASSOCIATED(helium%nnp_sr_cut)) THEN
     725      275860 :          DO i = 1, helium%nnp%num_atoms
     726      220688 :             IF (i == ind_he) CYCLE
     727      662064 :             rvect(:) = helium%nnp%coord(:, i) - helium%nnp%coord(:, ind_he)
     728      165516 :             CALL helium_pbc(helium, rvect)
     729      165516 :             rsqr = rvect(1)*rvect(1) + rvect(2)*rvect(2) + rvect(3)*rvect(3)
     730      220688 :             IF (rsqr < helium%nnp_sr_cut(helium%nnp%ele_ind(i))) THEN
     731           0 :                energy = 0.3_dp + 1.0_dp/rsqr
     732           0 :                IF (PRESENT(force)) THEN
     733           0 :                   force = 0.0_dp
     734             :                END IF
     735           0 :                RETURN
     736             :             END IF
     737             :          END DO
     738             :       END IF
     739             : 
     740             :       ! reset flag if there's an extrapolation to report:
     741       55172 :       helium%nnp%output_expol = .FALSE.
     742             : 
     743             :       ! calc atomic contribution to energy and force
     744             : !NOTE corresponds to nnp_force line with parallelization:
     745             : !DO i = istart, istart + mecalc - 1
     746      275860 :       DO i = 1, helium%nnp%num_atoms
     747             : 
     748             :          !determine index of atom type
     749      220688 :          ind = helium%nnp%ele_ind(i)
     750             : 
     751             :          !reset input nodes and grads of ele(ind):
     752     4744792 :          helium%nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
     753      220688 :          nnp => helium%nnp ! work around wrong INTENT of nnp_calc_acsf
     754      220688 :          IF (PRESENT(force)) THEN
     755      137600 :             helium%nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
     756       25600 :             ALLOCATE (dsymdxyz(3, helium%nnp%arc(ind)%n_nodes(1), helium%nnp%num_atoms))
     757       19200 :             ALLOCATE (denergydsym(helium%nnp%arc(ind)%n_nodes(1)))
     758     2131200 :             dsymdxyz(:, :, :) = 0.0_dp
     759        6400 :             CALL nnp_calc_acsf(nnp, i, dsymdxyz)
     760             :          ELSE
     761      214288 :             CALL nnp_calc_acsf(nnp, i)
     762             :          END IF
     763             : 
     764             :          ! input nodes filled, perform prediction:
     765      441376 :          DO i_com = 1, helium%nnp%n_committee !loop over committee members
     766             :             ! Predict energy
     767      220688 :             CALL nnp_predict(helium%nnp%arc(ind), helium%nnp, i_com)
     768      220688 :             helium%nnp%atomic_energy(i, i_com) = helium%nnp%arc(ind)%layer(helium%nnp%n_layer)%node(1) ! + helium%nnp%atom_energies(ind)
     769             : 
     770             :             !Gradients
     771      441376 :             IF (PRESENT(force)) THEN
     772             : 
     773      137600 :                denergydsym(:) = 0.0_dp
     774             : 
     775        6400 :                CALL nnp_gradients(helium%nnp%arc(ind), helium%nnp, i_com, denergydsym)
     776      137600 :                DO j = 1, helium%nnp%arc(ind)%n_nodes(1)
     777      662400 :                   DO k = 1, helium%nnp%num_atoms
     778     2230400 :                      DO m = 1, 3
     779             :                         helium%nnp%myforce(m, k, i_com) = helium%nnp%myforce(m, k, i_com) &
     780     2099200 :                                                           - denergydsym(j)*dsymdxyz(m, j, k)
     781             :                      END DO
     782             :                   END DO
     783             :                END DO
     784             : 
     785             :             END IF
     786             :          END DO ! end loop over committee members
     787             : 
     788             :          !deallocate memory
     789      275860 :          IF (PRESENT(force)) THEN
     790        6400 :             DEALLOCATE (denergydsym)
     791        6400 :             DEALLOCATE (dsymdxyz)
     792             :          END IF
     793             : 
     794             :       END DO ! end loop over num_atoms
     795             : 
     796             :       ! calculate energy:
     797      331032 :       helium%nnp%committee_energy(:) = SUM(helium%nnp%atomic_energy, 1)
     798      110344 :       energy = SUM(helium%nnp%committee_energy)/REAL(helium%nnp%n_committee, dp)
     799       55172 :       helium%nnp%nnp_potential_energy = energy
     800             : 
     801       55172 :       IF (PRESENT(force)) THEN
     802             :          ! bring myforce to force array
     803        8000 :          DO j = 1, helium%nnp%num_atoms
     804       27200 :             DO k = 1, 3
     805       44800 :                helium%nnp%committee_forces(k, j, :) = helium%nnp%myforce(k, j, :)
     806             :             END DO
     807             :          END DO
     808       46400 :          helium%nnp%nnp_forces(:, :) = SUM(helium%nnp%committee_forces, DIM=3)/REAL(helium%nnp%n_committee, dp)
     809             :          ! project out helium force entry
     810        1600 :          ig = 1
     811        8000 :          DO j = 1, helium%nnp%num_atoms
     812        6400 :             IF (j == ind_he) CYCLE
     813       19200 :             DO k = 1, 3
     814       19200 :                force(3*(helium%nnp%sort(ig) - 1) + k) = helium%nnp%nnp_forces(k, j)
     815             :             END DO
     816        8000 :             ig = ig + 1
     817             :          END DO
     818             :       END IF
     819             : 
     820             :       ! print properties if requested
     821       55172 :       print_section => section_vals_get_subs_vals(helium%nnp%nnp_input, "PRINT")
     822       55172 :       CALL helium_nnp_print(helium%nnp, print_section, ind_he)
     823             : 
     824       55172 :       RETURN
     825             : 
     826       55172 :    END SUBROUTINE helium_intpot_nnp
     827             : 
     828             : ! ***************************************************************************
     829             : !> \brief Helium-helium pair interaction potential.
     830             : !> \param r ...
     831             : !> \return ...
     832             : ! **************************************************************************************************
     833        1300 :    ELEMENTAL FUNCTION helium_vij(r) RESULT(vij)
     834             : 
     835             :       REAL(kind=dp), INTENT(IN)                          :: r
     836             :       REAL(kind=dp)                                      :: vij
     837             : 
     838             :       REAL(kind=dp)                                      :: f, x, x2
     839             : 
     840        1300 :       x = angstrom*r/2.9673_dp
     841        1300 :       IF (x < 1.241314_dp) THEN
     842         351 :          x2 = 1.241314_dp/x - 1.0_dp
     843         351 :          f = EXP(-x2*x2)
     844             :       ELSE
     845             :          f = 1.0_dp
     846             :       END IF
     847        1300 :       x2 = 1.0_dp/(x*x)
     848             :       vij = 10.8_dp/kelvin*(544850.4_dp*EXP(-13.353384_dp*x) - f* &
     849        1300 :                             ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2)
     850        1300 :    END FUNCTION helium_vij
     851             : 
     852             : #if 0
     853             : 
     854             :    ! this block is currently turned off
     855             : 
     856             : ! ***************************************************************************
     857             : !> \brief Helium-helium pair interaction potential's derivative.
     858             : !> \param r ...
     859             : !> \return ...
     860             : ! **************************************************************************************************
     861             :    ELEMENTAL FUNCTION helium_d_vij(r) RESULT(dvij)
     862             : 
     863             :       REAL(kind=dp), INTENT(IN)                          :: r
     864             :       REAL(kind=dp)                                      :: dvij
     865             : 
     866             :       REAL(kind=dp)                                      :: f, fp, x, x2, y
     867             : 
     868             :       x = angstrom*r/2.9673_dp
     869             :       x = r/2.9673_dp
     870             :       x2 = 1.0_dp/(x*x)
     871             :       IF (x < 1.241314_dp) THEN
     872             :          y = 1.241314_dp/x - 1.0_dp
     873             :          f = EXP(-y*y)
     874             :          fp = 2.0_dp*1.241314_dp*f*y* &
     875             :               ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2*x2
     876             :       ELSE
     877             :          f = 1.0_dp
     878             :          fp = 0.0_dp
     879             :       END IF
     880             : 
     881             :       dvij = angstrom*(10.8_dp/2.9673_dp)*( &
     882             :              (-13.353384_dp*544850.4_dp)*EXP(-13.353384_dp*x) - fp + &
     883             :              f*(((10.0_dp*0.1781_dp)*x2 + (8.0_dp*0.4253785_dp))*x2 + (6.0_dp*1.3732412_dp))* &
     884             :              x2*x2*x2/x)/(r*kelvin)
     885             :    END FUNCTION helium_d_vij
     886             : 
     887             : ! **************************************************************************************************
     888             : !> \brief ...
     889             : !> \param helium ...
     890             : !> \param n ...
     891             : !> \param i ...
     892             : !> \return ...
     893             : ! **************************************************************************************************
     894             :    FUNCTION helium_atom_action(helium, n, i) RESULT(res)
     895             : 
     896             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     897             :       INTEGER, INTENT(IN)                                :: n, i
     898             :       REAL(KIND=dp)                                      :: res
     899             : 
     900             :       INTEGER                                            :: c, j
     901             :       REAL(KIND=dp)                                      :: r(3), rp(3), s, t
     902             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     903             : 
     904             :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
     905             :       s = 0.0_dp
     906             :       t = 0.0_dp
     907             :       IF (n < helium%beads) THEN
     908             :          DO c = 1, 3
     909             :             r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
     910             :          END DO
     911             :          CALL helium_pbc(helium, r)
     912             :          t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     913             :          DO j = 1, i - 1
     914             :             DO c = 1, 3
     915             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     916             :                rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
     917             :             END DO
     918             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     919             :          END DO
     920             :          DO j = i + 1, helium%atoms
     921             :             DO c = 1, 3
     922             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     923             :                rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
     924             :             END DO
     925             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     926             :          END DO
     927             :       ELSE
     928             :          DO c = 1, 3
     929             :             r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
     930             :          END DO
     931             :          CALL helium_pbc(helium, r)
     932             :          t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     933             :          DO j = 1, i - 1
     934             :             DO c = 1, 3
     935             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     936             :                rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
     937             :             END DO
     938             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     939             :          END DO
     940             :          DO j = i + 1, helium%atoms
     941             :             DO c = 1, 3
     942             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     943             :                rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
     944             :             END DO
     945             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     946             :          END DO
     947             :       END IF
     948             :       t = t/(2.0_dp*helium%tau*helium%hb2m)
     949             :       s = s*0.5_dp
     950             :       res = s + t
     951             :       DEALLOCATE (work3)
     952             : 
     953             :    END FUNCTION helium_atom_action
     954             : 
     955             : ! **************************************************************************************************
     956             : !> \brief ...
     957             : !> \param helium ...
     958             : !> \param n ...
     959             : !> \return ...
     960             : ! **************************************************************************************************
     961             :    FUNCTION helium_link_action(helium, n) RESULT(res)
     962             : 
     963             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     964             :       INTEGER, INTENT(IN)                                :: n
     965             :       REAL(KIND=dp)                                      :: res
     966             : 
     967             :       INTEGER                                            :: c, i, j
     968             :       REAL(KIND=dp)                                      :: r(3), rp(3), s, t
     969             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     970             : 
     971             :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
     972             :       s = 0.0_dp
     973             :       t = 0.0_dp
     974             :       IF (n < helium%beads) THEN
     975             :          DO i = 1, helium%atoms
     976             :             DO c = 1, 3
     977             :                r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
     978             :             END DO
     979             :             CALL helium_pbc(helium, r)
     980             :             t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     981             :             DO j = 1, i - 1
     982             :                DO c = 1, 3
     983             :                   r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     984             :                   rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
     985             :                END DO
     986             :                s = s + helium_eval_expansion(helium, r, rp, work3)
     987             :             END DO
     988             :          END DO
     989             :       ELSE
     990             :          DO i = 1, helium%atoms
     991             :             DO c = 1, 3
     992             :                r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
     993             :             END DO
     994             :             CALL helium_pbc(helium, r)
     995             :             t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     996             :             DO j = 1, i - 1
     997             :                DO c = 1, 3
     998             :                   r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     999             :                   rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
    1000             :                END DO
    1001             :                s = s + helium_eval_expansion(helium, r, rp, work3)
    1002             :             END DO
    1003             :          END DO
    1004             :       END IF
    1005             :       t = t/(2.0_dp*helium%tau*helium%hb2m)
    1006             :       res = s + t
    1007             :       DEALLOCATE (work3)
    1008             : 
    1009             :    END FUNCTION helium_link_action
    1010             : 
    1011             : ! **************************************************************************************************
    1012             : !> \brief ...
    1013             : !> \param helium ...
    1014             : !> \return ...
    1015             : ! **************************************************************************************************
    1016             :    FUNCTION helium_total_action(helium) RESULT(res)
    1017             : 
    1018             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1019             :       REAL(KIND=dp)                                      :: res
    1020             : 
    1021             :       INTEGER                                            :: i
    1022             :       REAL(KIND=dp)                                      :: s
    1023             : 
    1024             :       s = 0.0_dp
    1025             :       DO i = 1, helium%beads
    1026             :          s = s + helium_link_action(helium, i)
    1027             :       END DO
    1028             :       res = s
    1029             : 
    1030             :    END FUNCTION helium_total_action
    1031             : 
    1032             : ! **************************************************************************************************
    1033             : !> \brief ...
    1034             : !> \param helium ...
    1035             : !> \param part ...
    1036             : !> \param ref_bead ...
    1037             : !> \param delta_bead ...
    1038             : !> \param d ...
    1039             : ! **************************************************************************************************
    1040             :    SUBROUTINE helium_delta_pos(helium, part, ref_bead, delta_bead, d)
    1041             : 
    1042             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
    1043             :       INTEGER, INTENT(IN)                                :: part, ref_bead, delta_bead
    1044             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: d
    1045             : 
    1046             :       INTEGER                                            :: b, bead, db, nbead, np, p
    1047             :       REAL(KIND=dp), DIMENSION(3)                        :: r
    1048             : 
    1049             :       b = helium%beads
    1050             : 
    1051             :       d(:) = 0.0_dp
    1052             :       IF (delta_bead > 0) THEN
    1053             :          bead = ref_bead
    1054             :          p = part
    1055             :          db = delta_bead
    1056             :          DO
    1057             :             IF (db < 1) EXIT
    1058             :             nbead = bead + 1
    1059             :             np = p
    1060             :             IF (nbead > b) THEN
    1061             :                nbead = nbead - b
    1062             :                np = helium%permutation(np)
    1063             :             END IF
    1064             :             r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
    1065             :             CALL helium_pbc(helium, r)
    1066             :             d(:) = d(:) + r(:)
    1067             :             bead = nbead
    1068             :             p = np
    1069             :             db = db - 1
    1070             :          END DO
    1071             :       ELSEIF (delta_bead < 0) THEN
    1072             :          bead = ref_bead
    1073             :          p = part
    1074             :          db = delta_bead
    1075             :          DO
    1076             :             IF (db >= 0) EXIT
    1077             :             nbead = bead - 1
    1078             :             np = p
    1079             :             IF (nbead < 1) THEN
    1080             :                nbead = nbead + b
    1081             :                np = helium%iperm(np)
    1082             :             END IF
    1083             :             r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
    1084             :             CALL helium_pbc(helium, r)
    1085             :             d(:) = d(:) + r(:)
    1086             :             bead = nbead
    1087             :             p = np
    1088             :             db = db + 1
    1089             :          END DO
    1090             :       END IF
    1091             :    END SUBROUTINE helium_delta_pos
    1092             : 
    1093             : #endif
    1094             : 
    1095             : END MODULE helium_interactions

Generated by: LCOV version 1.15