LCOV - code coverage report
Current view: top level - src/motion - pint_pile.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.5 % 57 55
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief  Methods to apply a simple Lagevin thermostat to PI runs.
      10              : !>         v_new = c1*vold + SQRT(kT/m)*c2*random
      11              : !> \author Felix Uhl
      12              : !> \par History
      13              : !>      10.2014 created [Felix Uhl]
      14              : ! **************************************************************************************************
      15              : MODULE pint_pile
      16              :    USE input_constants,                 ONLY: propagator_cmd
      17              :    USE input_section_types,             ONLY: section_vals_get,&
      18              :                                               section_vals_get_subs_vals,&
      19              :                                               section_vals_type,&
      20              :                                               section_vals_val_get
      21              :    USE kinds,                           ONLY: dp
      22              :    USE parallel_rng_types,              ONLY: GAUSSIAN,&
      23              :                                               rng_record_length,&
      24              :                                               rng_stream_type,&
      25              :                                               rng_stream_type_from_record
      26              :    USE pint_types,                      ONLY: normalmode_env_type,&
      27              :                                               pile_therm_type,&
      28              :                                               pint_env_type
      29              : #include "../base/base_uses.f90"
      30              : 
      31              :    IMPLICIT NONE
      32              : 
      33              :    PRIVATE
      34              : 
      35              :    PUBLIC :: pint_pile_step, &
      36              :              pint_pile_init, &
      37              :              pint_pile_release, &
      38              :              pint_calc_pile_energy
      39              : 
      40              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_pile'
      41              : 
      42              : CONTAINS
      43              : 
      44              : ! ***************************************************************************
      45              : !> \brief initializes the data for a pile run
      46              : !> \param pile_therm ...
      47              : !> \param pint_env ...
      48              : !> \param normalmode_env ...
      49              : !> \param section ...
      50              : !> \author Felix Uhl
      51              : ! **************************************************************************************************
      52          250 :    SUBROUTINE pint_pile_init(pile_therm, pint_env, normalmode_env, section)
      53              :       TYPE(pile_therm_type), INTENT(OUT)                 :: pile_therm
      54              :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
      55              :       TYPE(normalmode_env_type), POINTER                 :: normalmode_env
      56              :       TYPE(section_vals_type), POINTER                   :: section
      57              : 
      58              :       CHARACTER(LEN=rng_record_length)                   :: rng_record
      59              :       INTEGER                                            :: i, i_propagator, j, p
      60              :       LOGICAL                                            :: explicit
      61              :       REAL(KIND=dp)                                      :: dti2, ex
      62              :       REAL(KIND=dp), DIMENSION(3, 2)                     :: initial_seed
      63              :       TYPE(section_vals_type), POINTER                   :: pint_section, rng_section
      64              : 
      65           10 :       pint_env%e_pile = 0.0_dp
      66              :       pile_therm%thermostat_energy = 0.0_dp
      67              :       !Get input parameter
      68           10 :       CALL section_vals_val_get(section, "TAU", r_val=pile_therm%tau)
      69           10 :       CALL section_vals_val_get(section, "LAMBDA", r_val=pile_therm%lamb)
      70           10 :       CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=pile_therm%thermostat_energy)
      71           10 :       pint_section => section_vals_get_subs_vals(pint_env%input, "MOTION%PINT")
      72           10 :       CALL section_vals_val_get(pint_section, "PROPAGATOR", i_val=i_propagator)
      73              : 
      74           10 :       IF (i_propagator == propagator_cmd) THEN
      75            2 :          pile_therm%tau = 0.0_dp
      76              :       END IF
      77              : 
      78           10 :       p = pint_env%p
      79           10 :       dti2 = 0.5_dp*pint_env%dt
      80           30 :       ALLOCATE (pile_therm%c1(p))
      81           20 :       ALLOCATE (pile_therm%c2(p))
      82           20 :       ALLOCATE (pile_therm%g_fric(p))
      83           40 :       ALLOCATE (pile_therm%massfact(p, pint_env%ndim))
      84              :       !Initialize everything
      85              :       ! If tau is negative or zero the thermostat does not act on the centroid
      86              :       ! (TRPMD)
      87           10 :       IF (pile_therm%tau <= 0.0_dp) THEN
      88            2 :          pile_therm%g_fric(1) = 0.0_dp
      89              :       ELSE
      90            8 :          pile_therm%g_fric(1) = 1.0_dp/pile_therm%tau
      91              :       END IF
      92          132 :       DO i = 2, p
      93          132 :          pile_therm%g_fric(i) = 2.0_dp*pile_therm%lamb*SQRT(normalmode_env%lambda(i))
      94              :       END DO
      95          142 :       DO i = 1, p
      96          132 :          ex = -dti2*pile_therm%g_fric(i)
      97          132 :          pile_therm%c1(i) = EXP(ex)
      98          132 :          ex = pile_therm%c1(i)*pile_therm%c1(i)
      99          142 :          pile_therm%c2(i) = SQRT(1.0_dp - ex)
     100              :       END DO
     101         9298 :       DO j = 1, pint_env%ndim
     102        47278 :          DO i = 1, pint_env%p
     103        47268 :             pile_therm%massfact(i, j) = SQRT(pint_env%kT/pint_env%mass_fict(i, j))
     104              :          END DO
     105              :       END DO
     106              : 
     107              :       !prepare Random number generator
     108           10 :       NULLIFY (rng_section)
     109              :       rng_section => section_vals_get_subs_vals(section, &
     110           10 :                                                 subsection_name="RNG_INIT")
     111           10 :       CALL section_vals_get(rng_section, explicit=explicit)
     112           10 :       IF (explicit) THEN
     113              :          CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
     114            0 :                                    i_rep_val=1, c_val=rng_record)
     115              : 
     116            0 :          pile_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
     117              :       ELSE
     118           90 :          initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp)
     119              :          pile_therm%gaussian_rng_stream = rng_stream_type( &
     120              :                                           name="pile_rng_gaussian", distribution_type=GAUSSIAN, &
     121              :                                           extended_precision=.TRUE., &
     122           10 :                                           seed=initial_seed)
     123              :       END IF
     124              : 
     125           20 :    END SUBROUTINE pint_pile_init
     126              : 
     127              : ! **************************************************************************************************
     128              : !> \brief ...
     129              : !> \param vold ...
     130              : !> \param vnew ...
     131              : !> \param p ...
     132              : !> \param ndim ...
     133              : !> \param first_mode ...
     134              : !> \param masses ...
     135              : !> \param pile_therm ...
     136              : ! **************************************************************************************************
     137          224 :    SUBROUTINE pint_pile_step(vold, vnew, p, ndim, first_mode, masses, pile_therm)
     138              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vold, vnew
     139              :       INTEGER, INTENT(IN)                                :: p, ndim, first_mode
     140              :       REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: masses
     141              :       TYPE(pile_therm_type), POINTER                     :: pile_therm
     142              : 
     143              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pint_pile_step'
     144              : 
     145              :       INTEGER                                            :: handle, ibead, idim
     146              :       REAL(KIND=dp)                                      :: delta_ekin
     147              : 
     148          224 :       CALL timeset(routineN, handle)
     149          224 :       delta_ekin = 0.0_dp
     150        94220 :       DO idim = 1, ndim
     151       484028 :          DO ibead = first_mode, p
     152              :             vnew(ibead, idim) = pile_therm%c1(ibead)*vold(ibead, idim) + &
     153              :                                 pile_therm%massfact(ibead, idim)*pile_therm%c2(ibead)* &
     154       389808 :                                 pile_therm%gaussian_rng_stream%next()
     155              :             delta_ekin = delta_ekin + masses(ibead, idim)*( &
     156              :                          vnew(ibead, idim)*vnew(ibead, idim) - &
     157       483804 :                          vold(ibead, idim)*vold(ibead, idim))
     158              :          END DO
     159              :       END DO
     160          224 :       pile_therm%thermostat_energy = pile_therm%thermostat_energy - 0.5_dp*delta_ekin
     161              : 
     162          224 :       CALL timestop(handle)
     163          224 :    END SUBROUTINE pint_pile_step
     164              : 
     165              : ! ***************************************************************************
     166              : !> \brief releases the pile environment
     167              : !> \param pile_therm pile data to be released
     168              : !> \author Felix Uhl
     169              : ! **************************************************************************************************
     170           10 :    SUBROUTINE pint_pile_release(pile_therm)
     171              : 
     172              :       TYPE(pile_therm_type), INTENT(INOUT)               :: pile_therm
     173              : 
     174           10 :       DEALLOCATE (pile_therm%c1)
     175           10 :       DEALLOCATE (pile_therm%c2)
     176           10 :       DEALLOCATE (pile_therm%g_fric)
     177           10 :       DEALLOCATE (pile_therm%massfact)
     178              : 
     179           10 :    END SUBROUTINE pint_pile_release
     180              : 
     181              : ! ***************************************************************************
     182              : !> \brief returns the pile kinetic energy contribution
     183              : !> \param pint_env ...
     184              : !> \author Felix Uhl
     185              : ! **************************************************************************************************
     186          122 :    SUBROUTINE pint_calc_pile_energy(pint_env)
     187              :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     188              : 
     189          122 :       IF (ASSOCIATED(pint_env%pile_therm)) THEN
     190          122 :          pint_env%e_pile = pint_env%pile_therm%thermostat_energy
     191              :       END IF
     192              : 
     193          122 :    END SUBROUTINE pint_calc_pile_energy
     194              : END MODULE
        

Generated by: LCOV version 2.0-1