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