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 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         200 :    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           8 :       pint_env%e_pile = 0.0_dp
      66           8 :       pile_therm%thermostat_energy = 0.0_dp
      67             :       !Get input parameter
      68           8 :       CALL section_vals_val_get(section, "TAU", r_val=pile_therm%tau)
      69           8 :       CALL section_vals_val_get(section, "LAMBDA", r_val=pile_therm%lamb)
      70           8 :       CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=pile_therm%thermostat_energy)
      71           8 :       pint_section => section_vals_get_subs_vals(pint_env%input, "MOTION%PINT")
      72           8 :       CALL section_vals_val_get(pint_section, "PROPAGATOR", i_val=i_propagator)
      73             : 
      74           8 :       IF (i_propagator == propagator_cmd) THEN
      75           2 :          pile_therm%tau = 0.0_dp
      76             :       END IF
      77             : 
      78           8 :       p = pint_env%p
      79           8 :       dti2 = 0.5_dp*pint_env%dt
      80          24 :       ALLOCATE (pile_therm%c1(p))
      81          24 :       ALLOCATE (pile_therm%c2(p))
      82          24 :       ALLOCATE (pile_therm%g_fric(p))
      83          32 :       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           8 :       IF (pile_therm%tau <= 0.0_dp) THEN
      88           2 :          pile_therm%g_fric(1) = 0.0_dp
      89             :       ELSE
      90           6 :          pile_therm%g_fric(1) = 1.0_dp/pile_therm%tau
      91             :       END IF
      92         100 :       DO i = 2, p
      93         100 :          pile_therm%g_fric(i) = 2.0_dp*pile_therm%lamb*SQRT(normalmode_env%lambda(i))
      94             :       END DO
      95         108 :       DO i = 1, p
      96         100 :          ex = -dti2*pile_therm%g_fric(i)
      97         100 :          pile_therm%c1(i) = EXP(ex)
      98         100 :          ex = pile_therm%c1(i)*pile_therm%c1(i)
      99         108 :          pile_therm%c2(i) = SQRT(1.0_dp - ex)
     100             :       END DO
     101        9278 :       DO j = 1, pint_env%ndim
     102       46970 :          DO i = 1, pint_env%p
     103       46962 :             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           8 :       NULLIFY (rng_section)
     109             :       rng_section => section_vals_get_subs_vals(section, &
     110           8 :                                                 subsection_name="RNG_INIT")
     111           8 :       CALL section_vals_get(rng_section, explicit=explicit)
     112           8 :       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          72 :          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           8 :                                           seed=initial_seed)
     123             :       END IF
     124             : 
     125          16 :    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         204 :    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         204 :       CALL timeset(routineN, handle)
     149         204 :       delta_ekin = 0.0_dp
     150       94020 :       DO idim = 1, ndim
     151      480948 :          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      386928 :                                 pile_therm%gaussian_rng_stream%next()
     155             :             delta_ekin = delta_ekin + masses(ibead, idim)*( &
     156             :                          vnew(ibead, idim)*vnew(ibead, idim) - &
     157      480744 :                          vold(ibead, idim)*vold(ibead, idim))
     158             :          END DO
     159             :       END DO
     160         204 :       pile_therm%thermostat_energy = pile_therm%thermostat_energy - 0.5_dp*delta_ekin
     161             : 
     162         204 :       CALL timestop(handle)
     163         204 :    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           8 :    SUBROUTINE pint_pile_release(pile_therm)
     171             : 
     172             :       TYPE(pile_therm_type), INTENT(INOUT)               :: pile_therm
     173             : 
     174           8 :       DEALLOCATE (pile_therm%c1)
     175           8 :       DEALLOCATE (pile_therm%c2)
     176           8 :       DEALLOCATE (pile_therm%g_fric)
     177           8 :       DEALLOCATE (pile_therm%massfact)
     178             : 
     179           8 :    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         110 :    SUBROUTINE pint_calc_pile_energy(pint_env)
     187             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     188             : 
     189         110 :       IF (ASSOCIATED(pint_env%pile_therm)) THEN
     190         110 :          pint_env%e_pile = pint_env%pile_therm%thermostat_energy
     191             :       END IF
     192             : 
     193         110 :    END SUBROUTINE pint_calc_pile_energy
     194             : END MODULE