LCOV - code coverage report
Current view: top level - src/motion - helium_interactions.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:dc34ec9) Lines: 198 252 78.6 %
Date: 2023-03-24 20:09:49 Functions: 9 10 90.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief  Methods 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 helium_common,                   ONLY: helium_eval_chain,&
      16             :                                               helium_eval_expansion,&
      17             :                                               helium_pbc,&
      18             :                                               helium_spline
      19             :    USE helium_types,                    ONLY: e_id_interact,&
      20             :                                               e_id_kinetic,&
      21             :                                               e_id_potential,&
      22             :                                               e_id_thermo,&
      23             :                                               e_id_total,&
      24             :                                               e_id_virial,&
      25             :                                               helium_solvent_p_type,&
      26             :                                               helium_solvent_type
      27             :    USE input_constants,                 ONLY: helium_sampling_worm,&
      28             :                                               helium_solute_intpot_mwater,&
      29             :                                               helium_solute_intpot_none
      30             :    USE kinds,                           ONLY: dp
      31             :    USE physcon,                         ONLY: angstrom,&
      32             :                                               kelvin
      33             :    USE pint_types,                      ONLY: pint_env_type
      34             :    USE splines_types,                   ONLY: spline_data_type
      35             : #include "../base/base_uses.f90"
      36             : 
      37             :    IMPLICIT NONE
      38             : 
      39             :    PRIVATE
      40             : 
      41             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      42             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_interactions'
      43             : 
      44             :    PUBLIC :: helium_calc_energy
      45             :    PUBLIC :: helium_total_link_action
      46             :    PUBLIC :: helium_total_pair_action
      47             :    PUBLIC :: helium_total_inter_action
      48             :    PUBLIC :: helium_solute_e_f
      49             :    PUBLIC :: helium_bead_solute_e_f
      50             :    PUBLIC :: helium_intpot_scan
      51             :    PUBLIC :: helium_vij
      52             : 
      53             : CONTAINS
      54             : 
      55             : ! ***************************************************************************
      56             : !> \brief  Calculate the helium energy (including helium-solute interaction)
      57             : !> \param    helium     helium environment
      58             : !> \param    pint_env   path integral environment
      59             : !> \par History
      60             : !>         2009-06 moved I/O out from here [lwalewski]
      61             : !> \author hforbert
      62             : ! **************************************************************************************************
      63        7042 :    SUBROUTINE helium_calc_energy(helium, pint_env)
      64             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
      65             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
      66             : 
      67             :       INTEGER                                            :: b, bead, i, j, n
      68        7042 :       INTEGER, DIMENSION(:), POINTER                     :: perm
      69             :       LOGICAL                                            :: nperiodic
      70             :       REAL(KIND=dp)                                      :: a, cell_size, en, interac, kin, pot, &
      71             :                                                             rmax, rmin, vkin
      72        7042 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work2, work3
      73        7042 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
      74             :       REAL(KIND=dp), DIMENSION(3)                        :: r
      75        7042 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pos
      76             :       TYPE(spline_data_type), POINTER                    :: e0
      77             : 
      78        7042 :       pos => helium%pos
      79        7042 :       perm => helium%permutation
      80        7042 :       e0 => helium%e0
      81        7042 :       cell_size = 0.5_dp*helium%cell_size
      82        7042 :       nperiodic = .NOT. helium%periodic
      83        7042 :       n = helium%atoms
      84        7042 :       b = helium%beads
      85        7042 :       en = 0.0_dp
      86        7042 :       pot = 0.0_dp
      87        7042 :       rmin = 1.0e20_dp
      88        7042 :       rmax = 0.0_dp
      89             :       ALLOCATE (work(3, helium%beads + 1), &
      90             :                 work2(helium%beads + 1), &
      91       49294 :                 work3(SIZE(helium%uoffdiag, 1) + 1))
      92      214004 :       DO i = 1, n - 1
      93     3502716 :          DO j = i + 1, n
      94    71130344 :             DO bead = 1, b
      95   274655240 :                work(:, bead) = pos(:, i, bead) - pos(:, j, bead)
      96             :             END DO
      97    13154848 :             work(:, b + 1) = pos(:, perm(i), 1) - pos(:, perm(j), 1)
      98     3288712 :             en = en + helium_eval_chain(helium, work, b + 1, work2, work3, energy=.TRUE.)
      99    71337306 :             DO bead = 1, b
     100    67841632 :                a = work2(bead)
     101    67841632 :                IF (a < rmin) rmin = a
     102    67841632 :                IF (a > rmax) rmax = a
     103    71130344 :                IF ((a < cell_size) .OR. nperiodic) THEN
     104    63517637 :                   pot = pot + helium_spline(helium%vij, a)
     105             :                END IF
     106             :             END DO
     107             :          END DO
     108             :       END DO
     109        7042 :       DEALLOCATE (work, work2, work3)
     110        7042 :       pot = pot/b
     111        7042 :       en = en/b
     112             : 
     113             :       ! helium-solute interaction energy (all beads of all particles)
     114        7042 :       interac = 0.0_dp
     115        7042 :       IF (helium%solute_present) THEN
     116        3632 :          CALL helium_solute_e(pint_env, helium, interac)
     117             :       END IF
     118        7042 :       interac = interac/b
     119             : 
     120             : !TODO:
     121        7042 :       vkin = 0.0_dp
     122             : !   vkin = helium_virial_energy(helium)
     123             : 
     124        7042 :       kin = 0.0_dp
     125      221046 :       DO i = 1, n
     126      856016 :          r(:) = pos(:, i, b) - pos(:, perm(i), 1)
     127      214004 :          CALL helium_pbc(helium, r)
     128      214004 :          kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     129     4413186 :          DO bead = 2, b
     130    16768560 :             r(:) = pos(:, i, bead - 1) - pos(:, i, bead)
     131     4192140 :             CALL helium_pbc(helium, r)
     132     4406144 :             kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     133             :          END DO
     134             :       END DO
     135        7042 :       kin = 1.5_dp*n/helium%tau - 0.5*kin/(b*helium%tau**2*helium%hb2m)
     136             : 
     137             : ! TODO: move printing somewhere else ?
     138             : !   print *,"POT = ",(pot/n+helium%e_corr)*kelvin,"K"
     139             : !   print *,"INTERAC = ",interac*kelvin,"K"
     140             : !   print *,"RMIN= ",rmin*angstrom,"A"
     141             : !   print *,"RMAX= ",rmax*angstrom,"A"
     142             : !   print *,"EVIRIAL not valid!"
     143             : !   print *,"ETHERMO= ",((en+kin)/n+helium%e_corr)*kelvin,"K"
     144             : !   print *,"ECORR= ",helium%e_corr*kelvin,"K"
     145             : !!   kin = helium_total_action(helium)
     146             : !!   print *,"ACTION= ",kin
     147             : !   print *,"WINDING#= ",helium_calc_winding(helium)
     148             : 
     149        7042 :       helium%energy_inst(e_id_potential) = pot/n + helium%e_corr
     150        7042 :       helium%energy_inst(e_id_kinetic) = (en - pot + kin)/n
     151        7042 :       helium%energy_inst(e_id_interact) = interac
     152        7042 :       helium%energy_inst(e_id_thermo) = (en + kin)/n + helium%e_corr
     153        7042 :       helium%energy_inst(e_id_virial) = vkin ! 0.0_dp at the moment
     154        7042 :       helium%energy_inst(e_id_total) = helium%energy_inst(e_id_thermo)
     155             :       ! Once vkin is properly implemented, switch to:
     156             :       ! helium%energy_inst(e_id_total) = (en+vkin)/n+helium%e_corr
     157             : 
     158       14084 :    END SUBROUTINE helium_calc_energy
     159             : 
     160             : ! ***************************************************************************
     161             : !> \brief  Computes the total harmonic link action of the helium
     162             : !> \param helium ...
     163             : !> \return ...
     164             : !> \date   2016-05-03
     165             : !> \author Felix Uhl
     166             : ! **************************************************************************************************
     167         120 :    REAL(KIND=dp) FUNCTION helium_total_link_action(helium) RESULT(linkaction)
     168             : 
     169             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     170             : 
     171             :       INTEGER                                            :: iatom, ibead
     172         120 :       INTEGER, DIMENSION(:), POINTER                     :: perm
     173             :       REAL(KIND=dp), DIMENSION(3)                        :: r
     174             : 
     175         120 :       perm => helium%permutation
     176         120 :       linkaction = 0.0_dp
     177             : 
     178             :       ! Harmonic Link action
     179             :       ! (r(m-1) - r(m))**2/(4*lambda*tau)
     180        2172 :       DO ibead = 1, helium%beads - 1
     181       47586 :          DO iatom = 1, helium%atoms
     182      181656 :             r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, iatom, ibead + 1)
     183       45414 :             CALL helium_pbc(helium, r)
     184       47466 :             linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     185             :          END DO
     186             :       END DO
     187        2610 :       DO iatom = 1, helium%atoms
     188             :          ! choose last bead connection according to permutation table
     189        9960 :          r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, perm(iatom), 1)
     190        2490 :          CALL helium_pbc(helium, r)
     191        2610 :          linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
     192             :       END DO
     193         120 :       linkaction = linkaction/(2.0_dp*helium%tau*helium%hb2m)
     194             : 
     195         120 :    END FUNCTION helium_total_link_action
     196             : 
     197             : ! ***************************************************************************
     198             : !> \brief  Computes the total pair action of the helium
     199             : !> \param helium ...
     200             : !> \return ...
     201             : !> \date   2016-05-03
     202             : !> \author Felix Uhl
     203             : ! **************************************************************************************************
     204         120 :    REAL(KIND=dp) FUNCTION helium_total_pair_action(helium) RESULT(pairaction)
     205             : 
     206             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     207             : 
     208             :       INTEGER                                            :: iatom, ibead, jatom, opatom, patom
     209         120 :       INTEGER, DIMENSION(:), POINTER                     :: perm
     210         120 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     211             :       REAL(KIND=dp), DIMENSION(3)                        :: r, rp
     212             : 
     213         360 :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
     214         120 :       perm => helium%permutation
     215         120 :       pairaction = 0.0_dp
     216             : 
     217             :       ! He-He pair action
     218        2172 :       DO ibead = 1, helium%beads - 1
     219       45534 :          DO iatom = 1, helium%atoms - 1
     220      698706 :             DO jatom = iatom + 1, helium%atoms
     221     2613168 :                r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
     222     2613168 :                rp(:) = helium%pos(:, iatom, ibead + 1) - helium%pos(:, jatom, ibead + 1)
     223      696654 :                pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     224             :             END DO
     225             :          END DO
     226             :       END DO
     227             :       !Ensure right permutation for pair action of last and first beads.
     228        2490 :       DO iatom = 1, helium%atoms - 1
     229       37710 :          DO jatom = iatom + 1, helium%atoms
     230      140880 :             r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, jatom, helium%beads)
     231      140880 :             rp(:) = helium%pos(:, perm(iatom), 1) - helium%pos(:, perm(jatom), 1)
     232       37590 :             pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     233             :          END DO
     234             :       END DO
     235             : 
     236             :       ! correct for open worm configurations
     237         120 :       IF (.NOT. helium%worm_is_closed) THEN
     238             :          ! special treatment if double bead is first bead
     239           0 :          iatom = helium%worm_atom_idx
     240           0 :          IF (helium%worm_bead_idx == 1) THEN
     241             :             ! patom is the atom in front of the lone head bead
     242           0 :             patom = helium%iperm(iatom)
     243             :             ! go through all atoms
     244           0 :             DO jatom = 1, helium%atoms
     245           0 :                IF (jatom == helium%worm_atom_idx) CYCLE
     246           0 :                opatom = helium%iperm(jatom)
     247             :                ! subtract pair action for closed link
     248           0 :                r(:) = helium%pos(:, iatom, 1) - helium%pos(:, jatom, 1)
     249           0 :                rp(:) = helium%pos(:, patom, helium%beads) - helium%pos(:, opatom, helium%beads)
     250           0 :                pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
     251             :                ! and add corrected extra link
     252             :                ! rp stays the same
     253           0 :                r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, 1)
     254           0 :                pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     255             :             END DO
     256             :          ELSE
     257             :             ! bead stays constant
     258           0 :             ibead = helium%worm_bead_idx
     259             :             ! go through all atoms
     260           0 :             DO jatom = 1, helium%atoms
     261           0 :                IF (jatom == helium%worm_atom_idx) CYCLE
     262             :                ! subtract pair action for closed link
     263           0 :                r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
     264           0 :                rp(:) = helium%pos(:, iatom, ibead - 1) - helium%pos(:, jatom, ibead - 1)
     265           0 :                pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
     266             :                ! and add corrected extra link
     267             :                ! rp stays the same
     268           0 :                r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, ibead)
     269           0 :                pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
     270             :             END DO
     271             :          END IF
     272             :       END IF
     273         120 :       DEALLOCATE (work3)
     274             : 
     275         120 :    END FUNCTION helium_total_pair_action
     276             : 
     277             : ! ***************************************************************************
     278             : !> \brief  Computes the total interaction of the helium with the solute
     279             : !> \param pint_env ...
     280             : !> \param helium ...
     281             : !> \return ...
     282             : !> \date   2016-05-03
     283             : !> \author Felix Uhl
     284             : ! **************************************************************************************************
     285         120 :    REAL(KIND=dp) FUNCTION helium_total_inter_action(pint_env, helium) RESULT(interaction)
     286             : 
     287             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     288             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     289             : 
     290             :       INTEGER                                            :: iatom, ibead
     291             :       REAL(KIND=dp)                                      :: e
     292             : 
     293         120 :       interaction = 0.0_dp
     294             : 
     295             :       ! InterAction with solute
     296         120 :       IF (helium%solute_present) THEN
     297        1564 :          DO ibead = 1, helium%beads
     298       27068 :             DO iatom = 1, helium%atoms
     299             : 
     300             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     301       25504 :                                            iatom, ibead, helium%pos(:, iatom, ibead), e)
     302       26976 :                interaction = interaction + e
     303             :             END DO
     304             :          END DO
     305          92 :          IF (helium%sampling_method == helium_sampling_worm) THEN
     306          12 :             IF (.NOT. helium%worm_is_closed) THEN
     307             :                ! subtract half of tail bead interaction again
     308             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     309             :                                            helium%worm_atom_idx, helium%worm_bead_idx, &
     310           0 :                                            helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), e)
     311           0 :                interaction = interaction - 0.5_dp*e
     312             :                ! add half of head bead interaction
     313             :                CALL helium_bead_solute_e_f(pint_env, helium, &
     314             :                                            helium%worm_atom_idx, helium%worm_bead_idx, &
     315           0 :                                            helium%worm_xtra_bead, e)
     316           0 :                interaction = interaction + 0.5_dp*e
     317             :             END IF
     318             :          END IF
     319             :       END IF
     320             : 
     321         120 :       interaction = interaction*helium%tau
     322             : 
     323         120 :    END FUNCTION helium_total_inter_action
     324             : 
     325             : ! ***************************************************************************
     326             : !> \brief Calculate general helium-solute interaction energy (and forces)
     327             : !>        between one helium bead and the corresponding solute time slice.
     328             : !> \param pint_env           path integral environment
     329             : !> \param helium ...
     330             : !> \param helium_part_index  helium particle index
     331             : !> \param helium_slice_index helium time slice index
     332             : !> \param helium_r_opt       explicit helium bead coordinates (optional)
     333             : !> \param energy             calculated energy
     334             : !> \param force              calculated force (if requested)
     335             : !> \par History
     336             : !>         2019-09 Added multiple-time striding in imag. time [cschran]
     337             : !> \author Lukasz Walewski
     338             : ! **************************************************************************************************
     339     5509296 :    SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
     340             :                                      helium_slice_index, helium_r_opt, energy, force)
     341             : 
     342             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     343             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     344             :       INTEGER, INTENT(IN)                                :: helium_part_index, helium_slice_index
     345             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL  :: helium_r_opt
     346             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     347             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     348             :          OPTIONAL, POINTER                               :: force
     349             : 
     350             :       INTEGER                                            :: hbeads, hi, qi, stride
     351             :       REAL(KIND=dp), DIMENSION(3)                        :: helium_r
     352     5509296 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: my_force
     353             : 
     354     5509296 :       hbeads = helium%beads
     355             :       ! helium bead index that is invariant wrt the rotations
     356     5509296 :       hi = MOD(helium_slice_index - 1 + hbeads + helium%relrot, hbeads) + 1
     357             :       ! solute bead index that belongs to hi helium index
     358     5509296 :       qi = ((hi - 1)*pint_env%p)/hbeads + 1
     359             : 
     360             :       ! coordinates of the helium bead
     361     5509296 :       IF (PRESENT(helium_r_opt)) THEN
     362     1429598 :          helium_r(:) = helium_r_opt(:)
     363             :       ELSE
     364    16318792 :          helium_r(:) = helium%pos(:, helium_part_index, helium_slice_index)
     365             :       END IF
     366             : 
     367    11018592 :       SELECT CASE (helium%solute_interaction)
     368             : 
     369             :       CASE (helium_solute_intpot_mwater)
     370     5509296 :          IF (PRESENT(force)) THEN
     371    54652288 :             force(:, :) = 0.0_dp
     372     1171264 :             my_force => force(qi, :)
     373             :             CALL helium_intpot_model_water( &
     374             :                pint_env%x(qi, :), &
     375             :                helium, &
     376             :                helium_r, &
     377             :                energy, &
     378             :                my_force &
     379     1171264 :                )
     380             :          ELSE
     381             :             CALL helium_intpot_model_water( &
     382             :                pint_env%x(qi, :), &
     383             :                helium, &
     384             :                helium_r, &
     385             :                energy &
     386     4338032 :                )
     387             :          END IF
     388             : 
     389             :       CASE (helium_solute_intpot_none)
     390           0 :          energy = 0.0_dp
     391     5509296 :          IF (PRESENT(force)) THEN
     392           0 :             force(:, :) = 0.0_dp
     393             :          END IF
     394             : 
     395             :       CASE DEFAULT
     396             : 
     397             :       END SELECT
     398             : 
     399             :       ! Account for Imaginary time striding in forces:
     400     5509296 :       IF (PRESENT(force)) THEN
     401     1171264 :          IF (hbeads < pint_env%p) THEN
     402        3072 :             stride = pint_env%p/hbeads
     403      915456 :             force = force*REAL(stride, dp)
     404             :          END IF
     405             :       END IF
     406             : 
     407     5509296 :    END SUBROUTINE helium_bead_solute_e_f
     408             : 
     409             : ! ***************************************************************************
     410             : !> \brief Calculate total helium-solute interaction energy and forces.
     411             : !> \param   pint_env   path integral environment
     412             : !> \param helium ...
     413             : !> \param   energy     calculated interaction energy
     414             : !> \author Lukasz Walewski
     415             : ! **************************************************************************************************
     416        2642 :    SUBROUTINE helium_solute_e_f(pint_env, helium, energy)
     417             : 
     418             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     419             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     420             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     421             : 
     422             :       INTEGER                                            :: ia, ib, jb, jc
     423             :       REAL(KIND=dp)                                      :: my_energy
     424        2642 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: force
     425             : 
     426        2642 :       NULLIFY (force)
     427        2642 :       force => helium%force_inst
     428             : 
     429        2642 :       energy = 0.0_dp
     430      123044 :       force(:, :) = 0.0_dp
     431             : 
     432             :       ! calculate the total interaction energy and gradients between the
     433             :       ! solute and the helium, sum over all beads of all He particles
     434       75846 :       DO ia = 1, helium%atoms
     435     1247110 :          DO ib = 1, helium%beads
     436             :             CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, &
     437     1171264 :                                         energy=my_energy, force=helium%rtmp_p_ndim_2d)
     438     1171264 :             energy = energy + my_energy
     439     6015540 :             DO jb = 1, pint_env%p
     440    48881984 :                DO jc = 1, pint_env%ndim
     441    47710720 :                   force(jb, jc) = force(jb, jc) + helium%rtmp_p_ndim_2d(jb, jc)
     442             :                END DO
     443             :             END DO
     444             :          END DO
     445             :       END DO
     446             : 
     447        2642 :    END SUBROUTINE helium_solute_e_f
     448             : 
     449             : ! ***************************************************************************
     450             : !> \brief Calculate total helium-solute interaction energy.
     451             : !> \param   pint_env   path integral environment
     452             : !> \param helium ...
     453             : !> \param   energy     calculated interaction energy
     454             : !> \author Lukasz Walewski
     455             : ! **************************************************************************************************
     456        3632 :    SUBROUTINE helium_solute_e(pint_env, helium, energy)
     457             : 
     458             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     459             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     460             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     461             : 
     462             :       INTEGER                                            :: ia, ib
     463             :       REAL(KIND=dp)                                      :: my_energy
     464             : 
     465        3632 :       energy = 0.0_dp
     466             : 
     467      108516 :       DO ia = 1, helium%atoms
     468     1786660 :          DO ib = 1, helium%beads
     469     1678144 :             CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, energy=my_energy)
     470     1783028 :             energy = energy + my_energy
     471             :          END DO
     472             :       END DO
     473             : 
     474        3632 :    END SUBROUTINE helium_solute_e
     475             : 
     476             : ! ***************************************************************************
     477             : !> \brief  Scan the helium-solute interaction energy within the periodic cell
     478             : !> \param pint_env ...
     479             : !> \param helium_env ...
     480             : !> \date   2014-01-22
     481             : !> \par    History
     482             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     483             : !> \author Lukasz Walewski
     484             : ! **************************************************************************************************
     485           0 :    SUBROUTINE helium_intpot_scan(pint_env, helium_env)
     486             : 
     487             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     488             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     489             : 
     490             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_intpot_scan'
     491             : 
     492             :       INTEGER                                            :: handle, ic, ix, iy, iz, k, nbin
     493             :       LOGICAL                                            :: wrapped
     494             :       REAL(KIND=dp)                                      :: delr, my_en, ox, oy, oz
     495             :       REAL(kind=dp), DIMENSION(3)                        :: pbc1, pbc2, pos
     496             : 
     497           0 :       CALL timeset(routineN, handle)
     498             : 
     499             :       ! Perform scan only on ionode, since this is only used to output the intpot
     500           0 :       IF (pint_env%logger%para_env%is_source()) THEN
     501             :          ! Assume ionode always to have at least one helium_env
     502           0 :          k = 1
     503           0 :          helium_env(k)%helium%rho_inst(1, :, :, :) = 0.0_dp
     504           0 :          nbin = helium_env(k)%helium%rho_nbin
     505           0 :          delr = helium_env(k)%helium%rho_delr
     506           0 :          helium_env(k)%helium%center(:) = 0.0_dp
     507           0 :          ox = helium_env(k)%helium%center(1) - helium_env(k)%helium%rho_maxr/2.0_dp
     508           0 :          oy = helium_env(k)%helium%center(2) - helium_env(k)%helium%rho_maxr/2.0_dp
     509           0 :          oz = helium_env(k)%helium%center(3) - helium_env(k)%helium%rho_maxr/2.0_dp
     510             : 
     511           0 :          DO ix = 1, nbin
     512           0 :             DO iy = 1, nbin
     513           0 :                DO iz = 1, nbin
     514             : 
     515             :                   ! put the probe in the center of the current voxel
     516           0 :                   pos(:) = (/ox + (ix - 0.5_dp)*delr, oy + (iy - 0.5_dp)*delr, oz + (iz - 0.5_dp)*delr/)
     517             : 
     518             :                   ! calc interaction energy for the current probe position
     519           0 :                   helium_env(k)%helium%pos(:, 1, 1) = pos(:)
     520           0 :                   CALL helium_bead_solute_e_f(pint_env, helium_env(k)%helium, 1, 1, energy=my_en)
     521             : 
     522             :                   ! check if the probe fits within the unit cell
     523           0 :                   pbc1(:) = pos(:) - helium_env(k)%helium%center
     524           0 :                   pbc2(:) = pbc1(:)
     525           0 :                   CALL helium_pbc(helium_env(k)%helium, pbc2)
     526           0 :                   wrapped = .FALSE.
     527           0 :                   DO ic = 1, 3
     528           0 :                      IF (ABS(pbc1(ic) - pbc2(ic)) .GT. 10.0_dp*EPSILON(0.0_dp)) THEN
     529           0 :                         wrapped = .TRUE.
     530             :                      END IF
     531             :                   END DO
     532             : 
     533             :                   ! set the interaction energy value
     534           0 :                   IF (wrapped) THEN
     535           0 :                      helium_env(k)%helium%rho_inst(1, ix, iy, iz) = 0.0_dp
     536             :                   ELSE
     537           0 :                      helium_env(k)%helium%rho_inst(1, ix, iy, iz) = my_en
     538             :                   END IF
     539             : 
     540             :                END DO
     541             :             END DO
     542             :          END DO
     543             :       END IF
     544             : 
     545           0 :       CALL timestop(handle)
     546           0 :    END SUBROUTINE helium_intpot_scan
     547             : 
     548             : ! ***************************************************************************
     549             : !> \brief Calculate model helium-solute interaction energy and forces
     550             : !>        between one helium bead and the corresponding solute time
     551             : !>        slice asuming water solute.
     552             : !> \param solute_x  solute positions ARR(3*NATOMS)
     553             : !>        to global atom indices
     554             : !> \param helium    only needed for helium_pbc call at the moment
     555             : !> \param helium_x  helium bead position ARR(3)
     556             : !> \param energy    calculated interaction energy
     557             : !> \param force ...
     558             : !> \author Felix Uhl
     559             : ! **************************************************************************************************
     560     5509296 :    SUBROUTINE helium_intpot_model_water(solute_x, helium, helium_x, energy, force)
     561             : 
     562             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: solute_x
     563             :       TYPE(helium_solvent_type), INTENT(IN)              :: helium
     564             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: helium_x
     565             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     566             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
     567             :          OPTIONAL, POINTER                               :: force
     568             : 
     569             :       INTEGER                                            :: i, ig
     570             :       REAL(KIND=dp)                                      :: d, d2, dd, ep, eps, s1, s2, sig
     571             :       REAL(KIND=dp), DIMENSION(3)                        :: dr, solute_r
     572             : 
     573     5509296 :       energy = 0.0_dp
     574     5509296 :       IF (PRESENT(force)) THEN
     575    11712640 :          force(:) = 0.0_dp
     576             :       END IF
     577             : 
     578     5509296 :       sig = 2.69_dp ! 1.4 Angstrom
     579     5509296 :       eps = 60.61e-6_dp ! 19 K
     580     5509296 :       s1 = 0.0_dp
     581    22037184 :       DO i = 1, SIZE(helium%solute_element)
     582    22037184 :          IF (helium%solute_element(i) == "H ") THEN
     583    11018592 :             ig = i - 1
     584    11018592 :             solute_r(1) = solute_x(3*ig + 1)
     585    11018592 :             solute_r(2) = solute_x(3*ig + 2)
     586    11018592 :             solute_r(3) = solute_x(3*ig + 3)
     587    44074368 :             dr(:) = solute_r(:) - helium_x(:)
     588    11018592 :             CALL helium_pbc(helium, dr)
     589    11018592 :             d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
     590    11018592 :             d = SQRT(d2)
     591    11018592 :             dd = (sig/d)**6
     592    11018592 :             ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
     593    11018592 :             s1 = s1 + ep
     594    11018592 :             s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
     595    11018592 :             IF (PRESENT(force)) THEN
     596     2342528 :                force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
     597     2342528 :                force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
     598     2342528 :                force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
     599             :             END IF
     600             :          END IF
     601             :       END DO ! i = 1, num_hydrogen
     602     5509296 :       energy = energy + s1
     603             : 
     604     5509296 :       sig = 5.01_dp ! 2.6 Angstrom
     605     5509296 :       eps = 104.5e-6_dp ! 33 K
     606     5509296 :       s1 = 0.0_dp
     607    22037184 :       DO i = 1, SIZE(helium%solute_element)
     608    22037184 :          IF (helium%solute_element(i) == "O ") THEN
     609     5509296 :             ig = i - 1
     610     5509296 :             solute_r(1) = solute_x(3*ig + 1)
     611     5509296 :             solute_r(2) = solute_x(3*ig + 2)
     612     5509296 :             solute_r(3) = solute_x(3*ig + 3)
     613    22037184 :             dr(:) = solute_r(:) - helium_x(:)
     614     5509296 :             CALL helium_pbc(helium, dr)
     615     5509296 :             d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
     616     5509296 :             d = SQRT(d2)
     617     5509296 :             dd = (sig/d)**6
     618     5509296 :             ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
     619     5509296 :             s1 = s1 + ep
     620     5509296 :             s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
     621     5509296 :             IF (PRESENT(force)) THEN
     622     1171264 :                force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
     623     1171264 :                force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
     624     1171264 :                force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
     625             :             END IF
     626             :          END IF
     627             :       END DO ! i = 1, num_chlorine
     628     5509296 :       energy = energy + s1
     629             : 
     630     5509296 :    END SUBROUTINE helium_intpot_model_water
     631             : 
     632             : ! ***************************************************************************
     633             : !> \brief Helium-helium pair interaction potential.
     634             : !> \param r ...
     635             : !> \return ...
     636             : ! **************************************************************************************************
     637        1200 :    ELEMENTAL FUNCTION helium_vij(r) RESULT(vij)
     638             : 
     639             :       REAL(kind=dp), INTENT(IN)                          :: r
     640             :       REAL(kind=dp)                                      :: vij
     641             : 
     642             :       REAL(kind=dp)                                      :: f, x, x2
     643             : 
     644        1200 :       x = angstrom*r/2.9673_dp
     645        1200 :       IF (x < 1.241314_dp) THEN
     646         324 :          x2 = 1.241314_dp/x - 1.0_dp
     647         324 :          f = EXP(-x2*x2)
     648             :       ELSE
     649             :          f = 1.0_dp
     650             :       END IF
     651        1200 :       x2 = 1.0_dp/(x*x)
     652             :       vij = 10.8_dp/kelvin*(544850.4_dp*EXP(-13.353384_dp*x) - f* &
     653        1200 :                             ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2)
     654        1200 :    END FUNCTION helium_vij
     655             : 
     656             : #if 0
     657             : 
     658             :    ! this block is currently turned off
     659             : 
     660             : ! ***************************************************************************
     661             : !> \brief Helium-helium pair interaction potential's derivative.
     662             : !> \param r ...
     663             : !> \return ...
     664             : ! **************************************************************************************************
     665             :    ELEMENTAL FUNCTION helium_d_vij(r) RESULT(dvij)
     666             : 
     667             :       REAL(kind=dp), INTENT(IN)                          :: r
     668             :       REAL(kind=dp)                                      :: dvij
     669             : 
     670             :       REAL(kind=dp)                                      :: f, fp, x, x2, y
     671             : 
     672             :       x = angstrom*r/2.9673_dp
     673             :       x = r/2.9673_dp
     674             :       x2 = 1.0_dp/(x*x)
     675             :       IF (x < 1.241314_dp) THEN
     676             :          y = 1.241314_dp/x - 1.0_dp
     677             :          f = EXP(-y*y)
     678             :          fp = 2.0_dp*1.241314_dp*f*y* &
     679             :               ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2*x2
     680             :       ELSE
     681             :          f = 1.0_dp
     682             :          fp = 0.0_dp
     683             :       END IF
     684             : 
     685             :       dvij = angstrom*(10.8_dp/2.9673_dp)*( &
     686             :              (-13.353384_dp*544850.4_dp)*EXP(-13.353384_dp*x) - fp + &
     687             :              f*(((10.0_dp*0.1781_dp)*x2 + (8.0_dp*0.4253785_dp))*x2 + (6.0_dp*1.3732412_dp))* &
     688             :              x2*x2*x2/x)/(r*kelvin)
     689             :    END FUNCTION helium_d_vij
     690             : 
     691             : ! **************************************************************************************************
     692             : !> \brief ...
     693             : !> \param helium ...
     694             : !> \param n ...
     695             : !> \param i ...
     696             : !> \return ...
     697             : ! **************************************************************************************************
     698             :    FUNCTION helium_atom_action(helium, n, i) RESULT(res)
     699             : 
     700             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     701             :       INTEGER, INTENT(IN)                                :: n, i
     702             :       REAL(KIND=dp)                                      :: res
     703             : 
     704             :       INTEGER                                            :: c, j
     705             :       REAL(KIND=dp)                                      :: r(3), rp(3), s, t
     706             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     707             : 
     708             :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
     709             :       s = 0.0_dp
     710             :       t = 0.0_dp
     711             :       IF (n < helium%beads) THEN
     712             :          DO c = 1, 3
     713             :             r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
     714             :          END DO
     715             :          CALL helium_pbc(helium, r)
     716             :          t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     717             :          DO j = 1, i - 1
     718             :             DO c = 1, 3
     719             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     720             :                rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
     721             :             END DO
     722             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     723             :          END DO
     724             :          DO j = i + 1, helium%atoms
     725             :             DO c = 1, 3
     726             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     727             :                rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
     728             :             END DO
     729             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     730             :          END DO
     731             :       ELSE
     732             :          DO c = 1, 3
     733             :             r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
     734             :          END DO
     735             :          CALL helium_pbc(helium, r)
     736             :          t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     737             :          DO j = 1, i - 1
     738             :             DO c = 1, 3
     739             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     740             :                rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
     741             :             END DO
     742             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     743             :          END DO
     744             :          DO j = i + 1, helium%atoms
     745             :             DO c = 1, 3
     746             :                r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     747             :                rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
     748             :             END DO
     749             :             s = s + helium_eval_expansion(helium, r, rp, work3)
     750             :          END DO
     751             :       END IF
     752             :       t = t/(2.0_dp*helium%tau*helium%hb2m)
     753             :       s = s*0.5_dp
     754             :       res = s + t
     755             :       DEALLOCATE (work3)
     756             : 
     757             :    END FUNCTION helium_atom_action
     758             : 
     759             : ! **************************************************************************************************
     760             : !> \brief ...
     761             : !> \param helium ...
     762             : !> \param n ...
     763             : !> \return ...
     764             : ! **************************************************************************************************
     765             :    FUNCTION helium_link_action(helium, n) RESULT(res)
     766             : 
     767             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     768             :       INTEGER, INTENT(IN)                                :: n
     769             :       REAL(KIND=dp)                                      :: res
     770             : 
     771             :       INTEGER                                            :: c, i, j
     772             :       REAL(KIND=dp)                                      :: r(3), rp(3), s, t
     773             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work3
     774             : 
     775             :       ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
     776             :       s = 0.0_dp
     777             :       t = 0.0_dp
     778             :       IF (n < helium%beads) THEN
     779             :          DO i = 1, helium%atoms
     780             :             DO c = 1, 3
     781             :                r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
     782             :             END DO
     783             :             CALL helium_pbc(helium, r)
     784             :             t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     785             :             DO j = 1, i - 1
     786             :                DO c = 1, 3
     787             :                   r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     788             :                   rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
     789             :                END DO
     790             :                s = s + helium_eval_expansion(helium, r, rp, work3)
     791             :             END DO
     792             :          END DO
     793             :       ELSE
     794             :          DO i = 1, helium%atoms
     795             :             DO c = 1, 3
     796             :                r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
     797             :             END DO
     798             :             CALL helium_pbc(helium, r)
     799             :             t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
     800             :             DO j = 1, i - 1
     801             :                DO c = 1, 3
     802             :                   r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
     803             :                   rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
     804             :                END DO
     805             :                s = s + helium_eval_expansion(helium, r, rp, work3)
     806             :             END DO
     807             :          END DO
     808             :       END IF
     809             :       t = t/(2.0_dp*helium%tau*helium%hb2m)
     810             :       res = s + t
     811             :       DEALLOCATE (work3)
     812             : 
     813             :    END FUNCTION helium_link_action
     814             : 
     815             : ! **************************************************************************************************
     816             : !> \brief ...
     817             : !> \param helium ...
     818             : !> \return ...
     819             : ! **************************************************************************************************
     820             :    FUNCTION helium_total_action(helium) RESULT(res)
     821             : 
     822             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     823             :       REAL(KIND=dp)                                      :: res
     824             : 
     825             :       INTEGER                                            :: i
     826             :       REAL(KIND=dp)                                      :: s
     827             : 
     828             :       s = 0.0_dp
     829             :       DO i = 1, helium%beads
     830             :          s = s + helium_link_action(helium, i)
     831             :       END DO
     832             :       res = s
     833             : 
     834             :    END FUNCTION helium_total_action
     835             : 
     836             : ! **************************************************************************************************
     837             : !> \brief ...
     838             : !> \param helium ...
     839             : !> \param part ...
     840             : !> \param ref_bead ...
     841             : !> \param delta_bead ...
     842             : !> \param d ...
     843             : ! **************************************************************************************************
     844             :    SUBROUTINE helium_delta_pos(helium, part, ref_bead, delta_bead, d)
     845             : 
     846             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
     847             :       INTEGER, INTENT(IN)                                :: part, ref_bead, delta_bead
     848             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: d
     849             : 
     850             :       INTEGER                                            :: b, bead, db, nbead, np, p
     851             :       REAL(KIND=dp), DIMENSION(3)                        :: r
     852             : 
     853             :       b = helium%beads
     854             : 
     855             :       d(:) = 0.0_dp
     856             :       IF (delta_bead > 0) THEN
     857             :          bead = ref_bead
     858             :          p = part
     859             :          db = delta_bead
     860             :          DO
     861             :             IF (db < 1) EXIT
     862             :             nbead = bead + 1
     863             :             np = p
     864             :             IF (nbead > b) THEN
     865             :                nbead = nbead - b
     866             :                np = helium%permutation(np)
     867             :             END IF
     868             :             r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
     869             :             CALL helium_pbc(helium, r)
     870             :             d(:) = d(:) + r(:)
     871             :             bead = nbead
     872             :             p = np
     873             :             db = db - 1
     874             :          END DO
     875             :       ELSEIF (delta_bead < 0) THEN
     876             :          bead = ref_bead
     877             :          p = part
     878             :          db = delta_bead
     879             :          DO
     880             :             IF (db >= 0) EXIT
     881             :             nbead = bead - 1
     882             :             np = p
     883             :             IF (nbead < 1) THEN
     884             :                nbead = nbead + b
     885             :                np = helium%iperm(np)
     886             :             END IF
     887             :             r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
     888             :             CALL helium_pbc(helium, r)
     889             :             d(:) = d(:) + r(:)
     890             :             bead = nbead
     891             :             p = np
     892             :             db = db + 1
     893             :          END DO
     894             :       END IF
     895             :    END SUBROUTINE helium_delta_pos
     896             : 
     897             : #endif
     898             : 
     899             : END MODULE helium_interactions

Generated by: LCOV version 1.15