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 Module with utility to perform MD Nudged Elastic Band Calculation
10 : !> \note
11 : !> Numerical accuracy for parallel runs:
12 : !> Each replica starts the SCF run from the one optimized
13 : !> in a previous run. It may happen then energies and derivatives
14 : !> of a serial run and a parallel run could be slightly different
15 : !> 'cause of a different starting density matrix.
16 : !> Exact results are obtained using:
17 : !> EXTRAPOLATION USE_GUESS in QS section (Teo 09.2006)
18 : !> \author Teodoro Laino 10.2006
19 : ! **************************************************************************************************
20 : MODULE neb_md_utils
21 : USE cp_units, ONLY: cp_unit_from_cp2k
22 : USE global_types, ONLY: global_environment_type
23 : USE input_constants, ONLY: band_md_opt,&
24 : do_band_collective
25 : USE input_section_types, ONLY: section_vals_get,&
26 : section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: dp
30 : USE neb_types, ONLY: neb_type,&
31 : neb_var_type
32 : USE particle_types, ONLY: get_particle_pos_or_vel,&
33 : particle_type,&
34 : update_particle_pos_or_vel
35 : USE physcon, ONLY: kelvin,&
36 : massunit
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_md_utils'
42 :
43 : PUBLIC :: neb_initialize_velocity, &
44 : control_vels_a, &
45 : control_vels_b, &
46 : get_temperatures
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Initialize velocities of replica in an MD optimized algorithm within
52 : !> NEB
53 : !> \param vels ...
54 : !> \param neb_section ...
55 : !> \param particle_set ...
56 : !> \param i_rep ...
57 : !> \param iw ...
58 : !> \param globenv ...
59 : !> \param neb_env ...
60 : !> \par History
61 : !> 25.11.2010 Consider core-shell model (MK)
62 : !> \author Teodoro Laino 09.2006
63 : ! **************************************************************************************************
64 184 : SUBROUTINE neb_initialize_velocity(vels, neb_section, particle_set, i_rep, iw, &
65 : globenv, neb_env)
66 :
67 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vels
68 : TYPE(section_vals_type), POINTER :: neb_section
69 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
70 : INTEGER, INTENT(IN) :: i_rep, iw
71 : TYPE(global_environment_type), POINTER :: globenv
72 : TYPE(neb_type), POINTER :: neb_env
73 :
74 : INTEGER :: iatom, ivar, natom, nparticle, nvar
75 : REAL(KIND=dp) :: akin, mass, mass_tot, sc, temp, &
76 : temp_ext, tmp_r1
77 : REAL(KIND=dp), DIMENSION(3) :: v, vcom
78 : TYPE(section_vals_type), POINTER :: md_section
79 :
80 184 : IF (neb_env%opt_type == band_md_opt) THEN
81 70 : md_section => section_vals_get_subs_vals(neb_section, "OPTIMIZE_BAND%MD")
82 70 : CALL section_vals_val_get(md_section, "TEMPERATURE", r_val=temp_ext)
83 : ! Initialize velocity according to external temperature
84 70 : nparticle = SIZE(vels, 1)
85 70 : natom = SIZE(particle_set)
86 70 : mass_tot = 0.0_dp
87 70 : vcom(1:3) = 0.0_dp
88 70 : CALL globenv%gaussian_rng_stream%fill(vels(:, i_rep))
89 : ! Check always if BAND is working in Cartesian or in internal coordinates
90 : ! If working in cartesian coordinates let's get rid of the COM
91 : ! Compute also the total mass (both in Cartesian and internal)
92 70 : IF (neb_env%use_colvar) THEN
93 10 : nvar = nparticle
94 10 : mass_tot = REAL(nvar, KIND=dp)*massunit
95 : ELSE
96 1080 : DO iatom = 1, natom
97 1020 : mass = particle_set(iatom)%atomic_kind%mass
98 1020 : mass_tot = mass_tot + mass
99 1020 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, vels(:, i_rep))
100 4140 : vcom(1:3) = vcom(1:3) + mass*v(1:3)
101 : END DO
102 240 : vcom(1:3) = vcom(1:3)/mass_tot
103 : END IF
104 : ! Compute kinetic energy and temperature
105 70 : akin = 0.0_dp
106 70 : IF (neb_env%use_colvar) THEN
107 20 : DO ivar = 1, nvar
108 20 : akin = akin + 0.5_dp*massunit*vels(ivar, i_rep)*vels(ivar, i_rep)
109 : END DO
110 : ELSE
111 1080 : DO iatom = 1, natom
112 1020 : mass = particle_set(iatom)%atomic_kind%mass
113 4080 : v(1:3) = -vcom(1:3)
114 1020 : CALL update_particle_pos_or_vel(iatom, particle_set, v(1:3), vels(:, i_rep))
115 4140 : akin = akin + 0.5_dp*mass*DOT_PRODUCT(v(1:3), v(1:3))
116 : END DO
117 60 : nvar = 3*natom
118 : END IF
119 70 : temp = 2.0_dp*akin/REAL(nvar, KIND=dp)
120 : ! Scale velocities to get the correct initial temperature and
121 70 : sc = SQRT(temp_ext/temp)
122 3140 : vels(:, i_rep) = vels(:, i_rep)*sc
123 : ! Re-compute kinetic energya and temperature and velocity of COM
124 70 : akin = 0.0_dp
125 70 : vcom = 0.0_dp
126 70 : IF (neb_env%use_colvar) THEN
127 20 : DO ivar = 1, nvar
128 20 : akin = akin + 0.5_dp*massunit*vels(ivar, i_rep)*vels(ivar, i_rep)
129 : END DO
130 : ELSE
131 1080 : DO iatom = 1, natom
132 1020 : mass = particle_set(iatom)%atomic_kind%mass
133 1020 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, vels(:, i_rep))
134 4080 : vcom(1:3) = vcom(1:3) + mass*v(1:3)
135 4140 : akin = akin + 0.5_dp*mass*DOT_PRODUCT(v(1:3), v(1:3))
136 : END DO
137 : END IF
138 280 : vcom(1:3) = vcom(1:3)/mass_tot
139 : ! Dump information
140 70 : IF (iw > 0) THEN
141 5 : temp = 2.0_dp*akin/REAL(nvar, KIND=dp)
142 5 : tmp_r1 = cp_unit_from_cp2k(temp, "K")
143 : WRITE (iw, '(A,T61,F18.2,A2)') &
144 5 : ' NEB| Initial Temperature ', tmp_r1, " K"
145 : WRITE (iw, '(A,T61,F20.12)') &
146 5 : ' NEB| Centre of mass velocity in direction x:', vcom(1), &
147 5 : ' NEB| Centre of mass velocity in direction y:', vcom(2), &
148 10 : ' NEB| Centre of mass velocity in direction z:', vcom(3)
149 5 : WRITE (iw, '(T2,"NEB|",75("*"))')
150 : END IF
151 : ELSE
152 32490 : vels(:, i_rep) = 0.0_dp
153 : END IF
154 :
155 184 : END SUBROUTINE neb_initialize_velocity
156 :
157 : ! **************************************************************************************************
158 : !> \brief Control on velocities - I part
159 : !> \param vels ...
160 : !> \param particle_set ...
161 : !> \param tc_section ...
162 : !> \param vc_section ...
163 : !> \param output_unit ...
164 : !> \param istep ...
165 : !> \author Teodoro Laino 09.2006
166 : ! **************************************************************************************************
167 160 : SUBROUTINE control_vels_a(vels, particle_set, tc_section, vc_section, &
168 : output_unit, istep)
169 : TYPE(neb_var_type), POINTER :: vels
170 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
171 : TYPE(section_vals_type), POINTER :: tc_section, vc_section
172 : INTEGER, INTENT(IN) :: output_unit, istep
173 :
174 : INTEGER :: i, temp_tol_steps
175 : LOGICAL :: explicit
176 : REAL(KIND=dp) :: ext_temp, f_annealing, scale, temp_tol, &
177 : temploc, tmp_r1, tmp_r2
178 160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temperatures
179 :
180 : ! Temperature control
181 :
182 160 : CALL section_vals_get(tc_section, explicit=explicit)
183 160 : IF (explicit) THEN
184 60 : CALL section_vals_val_get(tc_section, "TEMP_TOL_STEPS", i_val=temp_tol_steps)
185 60 : CALL section_vals_val_get(tc_section, "TEMPERATURE", r_val=ext_temp)
186 60 : CALL section_vals_val_get(tc_section, "TEMP_TOL", r_val=temp_tol)
187 180 : ALLOCATE (temperatures(SIZE(vels%wrk, 2)))
188 : ! Computes temperatures
189 60 : CALL get_temperatures(vels, particle_set, temperatures, factor=1.0_dp)
190 : ! Possibly rescale
191 60 : IF (istep <= temp_tol_steps) THEN
192 260 : DO i = 2, SIZE(vels%wrk, 2) - 1
193 220 : temploc = temperatures(i)
194 260 : IF (ABS(temploc - ext_temp) > temp_tol) THEN
195 158 : IF (output_unit > 0) THEN
196 79 : tmp_r1 = cp_unit_from_cp2k(temploc, "K")
197 79 : tmp_r2 = cp_unit_from_cp2k(ext_temp, "K")
198 : WRITE (output_unit, '(T2,"NEB| Replica Nr.",I5,'// &
199 : '" - Velocity rescaled from: ",F12.6," to: ",F12.6,".")') &
200 79 : i, tmp_r1, tmp_r2
201 :
202 : END IF
203 158 : scale = SQRT(ext_temp/temploc)
204 8216 : vels%wrk(:, i) = scale*vels%wrk(:, i)
205 : END IF
206 : END DO
207 : END IF
208 120 : DEALLOCATE (temperatures)
209 : END IF
210 : ! Annealing
211 160 : CALL section_vals_get(vc_section, explicit=explicit)
212 160 : IF (explicit) THEN
213 160 : CALL section_vals_val_get(vc_section, "ANNEALING", r_val=f_annealing)
214 740 : DO i = 2, SIZE(vels%wrk, 2) - 1
215 27320 : vels%wrk(:, i) = f_annealing*vels%wrk(:, i)
216 : END DO
217 : END IF
218 160 : END SUBROUTINE control_vels_a
219 :
220 : ! **************************************************************************************************
221 : !> \brief Control on velocities - II part
222 : !> \param vels ...
223 : !> \param forces ...
224 : !> \param vc_section ...
225 : !> \author Teodoro Laino 09.2006
226 : ! **************************************************************************************************
227 160 : SUBROUTINE control_vels_b(vels, forces, vc_section)
228 : TYPE(neb_var_type), POINTER :: vels, forces
229 : TYPE(section_vals_type), POINTER :: vc_section
230 :
231 : INTEGER :: i
232 : LOGICAL :: explicit, lval
233 : REAL(KIND=dp) :: factor, norm
234 :
235 : ! Check the sign of V.dot.F
236 :
237 160 : CALL section_vals_get(vc_section, explicit=explicit)
238 160 : IF (explicit) THEN
239 160 : CALL section_vals_val_get(vc_section, "PROJ_VELOCITY_VERLET", l_val=lval)
240 160 : IF (lval) THEN
241 560 : DO i = 2, SIZE(vels%wrk, 2) - 1
242 18840 : norm = DOT_PRODUCT(forces%wrk(:, i), forces%wrk(:, i))
243 18840 : factor = DOT_PRODUCT(vels%wrk(:, i), forces%wrk(:, i))
244 560 : IF (factor > 0 .AND. (norm >= EPSILON(0.0_dp))) THEN
245 36608 : vels%wrk(:, i) = factor/norm*forces%wrk(:, i)
246 : ELSE
247 536 : vels%wrk(:, i) = 0.0_dp
248 : END IF
249 : END DO
250 : END IF
251 160 : CALL section_vals_val_get(vc_section, "SD_LIKE", l_val=lval)
252 160 : IF (lval) THEN
253 0 : DO i = 2, SIZE(vels%wrk, 2) - 1
254 0 : vels%wrk(:, i) = 0.0_dp
255 : END DO
256 : END IF
257 : END IF
258 160 : END SUBROUTINE control_vels_b
259 :
260 : ! **************************************************************************************************
261 : !> \brief Computes temperatures
262 : !> \param vels ...
263 : !> \param particle_set ...
264 : !> \param temperatures ...
265 : !> \param ekin ...
266 : !> \param factor ...
267 : !> \par History
268 : !> 24.11.2010 rewritten to include core-shell model (MK)
269 : !> \author Teodoro Laino 09.2006
270 : ! **************************************************************************************************
271 349 : SUBROUTINE get_temperatures(vels, particle_set, temperatures, ekin, factor)
272 :
273 : TYPE(neb_var_type), POINTER :: vels
274 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
275 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: temperatures
276 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: ekin
277 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: factor
278 :
279 : INTEGER :: i_rep, iatom, ivar, n_rep, natom, nvar
280 : REAL(KIND=dp) :: akin, mass, myfactor, temploc
281 : REAL(KIND=dp), DIMENSION(3) :: v
282 :
283 349 : myfactor = kelvin
284 349 : IF (PRESENT(factor)) myfactor = factor
285 2425 : IF (PRESENT(ekin)) ekin(:) = 0.0_dp
286 2536 : temperatures(:) = 0.0_dp
287 349 : nvar = SIZE(vels%wrk, 1)
288 349 : n_rep = SIZE(vels%wrk, 2)
289 349 : natom = SIZE(particle_set)
290 1838 : DO i_rep = 2, n_rep - 1
291 1489 : akin = 0.0_dp
292 1489 : IF (vels%in_use == do_band_collective) THEN
293 72 : DO ivar = 1, nvar
294 72 : akin = akin + 0.5_dp*massunit*vels%wrk(ivar, i_rep)*vels%wrk(ivar, i_rep)
295 : END DO
296 : ELSE
297 44832 : DO iatom = 1, natom
298 43379 : mass = particle_set(iatom)%atomic_kind%mass
299 43379 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, vels%wrk(:, i_rep))
300 174969 : akin = akin + 0.5_dp*mass*DOT_PRODUCT(v(1:3), v(1:3))
301 : END DO
302 1453 : nvar = 3*natom
303 : END IF
304 1489 : IF (PRESENT(ekin)) ekin(i_rep) = akin
305 1489 : temploc = 2.0_dp*akin/REAL(nvar, KIND=dp)
306 1838 : temperatures(i_rep) = temploc*myfactor
307 : END DO
308 :
309 349 : END SUBROUTINE get_temperatures
310 :
311 : END MODULE neb_md_utils
|
