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