1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Methods to apply a simple Lagevin thermostat to PI runs. 8!> v_new = c1*vold + SQRT(kT/m)*c2*random 9!> \author Felix Uhl 10!> \par History 11!> 10.2014 created [Felix Uhl] 12! ************************************************************************************************** 13MODULE pint_pile 14 USE input_section_types, ONLY: section_vals_get,& 15 section_vals_get_subs_vals,& 16 section_vals_type,& 17 section_vals_val_get 18 USE kinds, ONLY: dp 19 USE parallel_rng_types, ONLY: GAUSSIAN,& 20 create_rng_stream,& 21 delete_rng_stream,& 22 next_random_number,& 23 read_rng_stream,& 24 rng_record_length 25 USE pint_types, ONLY: normalmode_env_type,& 26 pile_therm_type,& 27 pint_env_type 28#include "../base/base_uses.f90" 29 30 IMPLICIT NONE 31 32 PRIVATE 33 34 PUBLIC :: pint_pile_step, & 35 pint_pile_init, & 36 pint_pile_release, & 37 pint_calc_pile_energy 38 39 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_pile' 40 41CONTAINS 42 43 ! *************************************************************************** 44 !> \brief initializes the data for a pile run 45 !> \author Felix Uhl 46 ! *************************************************************************** 47! ************************************************************************************************** 48!> \brief ... 49!> \param pile_therm ... 50!> \param pint_env ... 51!> \param normalmode_env ... 52!> \param section ... 53! ************************************************************************************************** 54 SUBROUTINE pint_pile_init(pile_therm, pint_env, normalmode_env, section) 55 TYPE(pile_therm_type), POINTER :: pile_therm 56 TYPE(pint_env_type), POINTER :: pint_env 57 TYPE(normalmode_env_type), POINTER :: normalmode_env 58 TYPE(section_vals_type), POINTER :: section 59 60 CHARACTER(len=*), PARAMETER :: routineN = 'pint_pile_init', routineP = moduleN//':'//routineN 61 62 CHARACTER(LEN=rng_record_length) :: rng_record 63 INTEGER :: i, j, p 64 LOGICAL :: explicit 65 REAL(KIND=dp) :: dti2, ex 66 REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed 67 TYPE(section_vals_type), POINTER :: rng_section 68 69 pint_env%e_pile = 0.0_dp 70 ALLOCATE (pile_therm) 71 pile_therm%ref_count = 1 72 pile_therm%thermostat_energy = 0.0_dp 73 !Get input parameter 74 CALL section_vals_val_get(section, "TAU", r_val=pile_therm%tau) 75 CALL section_vals_val_get(section, "LAMBDA", r_val=pile_therm%lamb) 76 77 p = pint_env%p 78 dti2 = 0.5_dp*pint_env%dt 79 ALLOCATE (pile_therm%c1(p)) 80 ALLOCATE (pile_therm%c2(p)) 81 ALLOCATE (pile_therm%g_fric(p)) 82 ALLOCATE (pile_therm%massfact(p, pint_env%ndim)) 83 !Initialize everything 84 ! If tau is negative or zero the thermostat does not act on the centroid 85 ! (TRPMD) 86 IF (pile_therm%tau <= 0.0_dp) THEN 87 pile_therm%g_fric(1) = 0.0_dp 88 ELSE 89 pile_therm%g_fric(1) = 1.0_dp/pile_therm%tau 90 END IF 91 DO i = 2, p 92 pile_therm%g_fric(i) = 2.0_dp*pile_therm%lamb*SQRT(normalmode_env%lambda(i)) 93 END DO 94 DO i = 1, p 95 ex = -dti2*pile_therm%g_fric(i) 96 pile_therm%c1(i) = EXP(ex) 97 ex = pile_therm%c1(i)*pile_therm%c1(i) 98 pile_therm%c2(i) = SQRT(1.0_dp - ex) 99 END DO 100 DO j = 1, pint_env%ndim 101 DO i = 1, pint_env%p 102 pile_therm%massfact(i, j) = SQRT(pint_env%kT/pint_env%mass_fict(i, j)) 103 END DO 104 END DO 105 106 !prepare Random number generator 107 NULLIFY (rng_section) 108 NULLIFY (pile_therm%gaussian_rng_stream) 109 rng_section => section_vals_get_subs_vals(section, & 110 subsection_name="RNG_INIT") 111 CALL section_vals_get(rng_section, explicit=explicit) 112 IF (explicit) THEN 113 CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", & 114 i_rep_val=1, c_val=rng_record) 115 CALL read_rng_stream(rng_stream=pile_therm%gaussian_rng_stream, & 116 rng_record=rng_record) 117 ELSE 118 initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp) 119 CALL create_rng_stream(rng_stream=pile_therm%gaussian_rng_stream, & 120 name="pile_rng_gaussian", distribution_type=GAUSSIAN, & 121 extended_precision=.TRUE., & 122 seed=initial_seed) 123 END IF 124 125 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 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', routineP = moduleN//':'//routineN 144 145 INTEGER :: handle, ibead, idim 146 REAL(KIND=dp) :: delta_ekin 147 148 CALL timeset(routineN, handle) 149 delta_ekin = 0.0_dp 150 DO idim = 1, ndim 151 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 next_random_number(pile_therm%gaussian_rng_stream) 155 delta_ekin = delta_ekin + masses(ibead, idim)*( & 156 vnew(ibead, idim)*vnew(ibead, idim) - & 157 vold(ibead, idim)*vold(ibead, idim)) 158 END DO 159 END DO 160 pile_therm%thermostat_energy = pile_therm%thermostat_energy - 0.5_dp*delta_ekin 161 162 CALL timestop(handle) 163 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! ************************************************************************************************** 171!> \brief ... 172!> \param pile_therm ... 173! ************************************************************************************************** 174 SUBROUTINE pint_pile_release(pile_therm) 175 176 TYPE(pile_therm_type), POINTER :: pile_therm 177 178 CHARACTER(len=*), PARAMETER :: routineN = 'pint_pile_release', & 179 routineP = moduleN//':'//routineN 180 181 IF (ASSOCIATED(pile_therm)) THEN 182 pile_therm%ref_count = pile_therm%ref_count - 1 183 IF (pile_therm%ref_count == 0) THEN 184 DEALLOCATE (pile_therm%c1) 185 DEALLOCATE (pile_therm%c2) 186 DEALLOCATE (pile_therm%g_fric) 187 DEALLOCATE (pile_therm%massfact) 188 CALL delete_rng_stream(pile_therm%gaussian_rng_stream) 189 DEALLOCATE (pile_therm) 190 END IF 191 END IF 192 NULLIFY (pile_therm) 193 194 RETURN 195 END SUBROUTINE pint_pile_release 196 197 ! *************************************************************************** 198 !> \brief returns the pile kinetic energy contribution 199 !> \author Felix Uhl 200 ! *************************************************************************** 201! ************************************************************************************************** 202!> \brief ... 203!> \param pint_env ... 204! ************************************************************************************************** 205 SUBROUTINE pint_calc_pile_energy(pint_env) 206 TYPE(pint_env_type), POINTER :: pint_env 207 208 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_pile_energy', & 209 routineP = moduleN//':'//routineN 210 211 IF (ASSOCIATED(pint_env%pile_therm)) THEN 212 pint_env%e_pile = pint_env%pile_therm%thermostat_energy 213 END IF 214 215 RETURN 216 217 END SUBROUTINE pint_calc_pile_energy 218END MODULE 219