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 Collection of utilities for setting-up and handle velocities in MD
10 : !> runs
11 : !> \author CJM
12 : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008
13 : !> reorganization of the original routines/modules
14 : !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
15 : !> (patch by Marcel Baer)
16 : ! **************************************************************************************************
17 : MODULE md_vel_utils
18 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind,&
21 : get_atomic_kind_set
22 : USE bibliography, ONLY: West2006,&
23 : cite_reference
24 : USE cell_types, ONLY: &
25 : cell_transform_input_cartesian, cell_type, use_perd_none, use_perd_x, use_perd_xy, &
26 : use_perd_xyz, use_perd_xz, use_perd_y, use_perd_yz, use_perd_z
27 : USE cp_linked_list_input, ONLY: cp_sll_val_next,&
28 : cp_sll_val_type
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
32 : cp_print_key_unit_nr
33 : USE cp_subsys_types, ONLY: cp_subsys_get,&
34 : cp_subsys_type
35 : USE cp_units, ONLY: cp_unit_from_cp2k
36 : USE extended_system_types, ONLY: npt_info_type
37 : USE force_env_methods, ONLY: force_env_calc_energy_force
38 : USE force_env_types, ONLY: force_env_get,&
39 : force_env_type
40 : USE force_env_utils, ONLY: force_env_rattle,&
41 : force_env_shake
42 : USE global_types, ONLY: global_environment_type
43 : USE input_constants, ONLY: &
44 : md_init_default, md_init_vib, npe_f_ensemble, npe_i_ensemble, &
45 : nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
46 : npt_ia_ensemble, reftraj_ensemble
47 : USE input_cp2k_binary_restarts, ONLY: read_binary_velocities
48 : USE input_restart_force_eval, ONLY: update_subsys
49 : USE input_section_types, ONLY: section_vals_get,&
50 : section_vals_get_subs_vals,&
51 : section_vals_list_get,&
52 : section_vals_type,&
53 : section_vals_val_get
54 : USE input_val_types, ONLY: val_get,&
55 : val_type
56 : USE kinds, ONLY: default_string_length,&
57 : dp
58 : USE mathconstants, ONLY: pi
59 : USE mathlib, ONLY: diamat_all
60 : USE md_ener_types, ONLY: md_ener_type
61 : USE md_environment_types, ONLY: get_md_env,&
62 : md_environment_type
63 : USE md_util, ONLY: read_vib_eigs_unformatted
64 : USE message_passing, ONLY: mp_para_env_type
65 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
66 : USE molecule_kind_types, ONLY: fixd_constraint_type,&
67 : get_molecule_kind,&
68 : get_molecule_kind_set,&
69 : molecule_kind_type
70 : USE parallel_rng_types, ONLY: UNIFORM,&
71 : rng_stream_type
72 : USE particle_list_types, ONLY: particle_list_type
73 : USE particle_types, ONLY: particle_type
74 : USE physcon, ONLY: kelvin
75 : USE shell_opt, ONLY: optimize_shell_core
76 : USE shell_potential_types, ONLY: shell_kind_type
77 : USE simpar_types, ONLY: simpar_type
78 : USE thermal_region_types, ONLY: thermal_region_type,&
79 : thermal_regions_type
80 : #include "../base/base_uses.f90"
81 :
82 : IMPLICIT NONE
83 :
84 : PRIVATE
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_vel_utils'
86 :
87 : PUBLIC :: temperature_control, &
88 : comvel_control, &
89 : angvel_control, &
90 : setup_velocities
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief compute center of mass position
96 : !> *** is only used by initialize_velocities below ***
97 : !> \param part ...
98 : !> \param is_fixed ...
99 : !> \param rcom ...
100 : !> \par History
101 : !> 2007-11-6: created
102 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
103 : ! **************************************************************************************************
104 108 : SUBROUTINE compute_rcom(part, is_fixed, rcom)
105 : TYPE(particle_type), DIMENSION(:), POINTER :: part
106 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
107 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: rcom
108 :
109 : INTEGER :: i
110 : REAL(KIND=dp) :: denom, mass
111 : TYPE(atomic_kind_type), POINTER :: atomic_kind
112 :
113 108 : rcom(:) = 0.0_dp
114 108 : denom = 0.0_dp
115 1269 : DO i = 1, SIZE(part)
116 1161 : atomic_kind => part(i)%atomic_kind
117 1161 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
118 1269 : SELECT CASE (is_fixed(i))
119 : CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
120 1161 : rcom(1) = rcom(1) + part(i)%r(1)*mass
121 1161 : rcom(2) = rcom(2) + part(i)%r(2)*mass
122 1161 : rcom(3) = rcom(3) + part(i)%r(3)*mass
123 1161 : denom = denom + mass
124 : END SELECT
125 : END DO
126 432 : rcom = rcom/denom
127 :
128 108 : END SUBROUTINE compute_rcom
129 :
130 : ! **************************************************************************************************
131 : !> \brief compute center of mass velocity
132 : !> *** is only used by initialize_velocities below ***
133 : !> \param part ...
134 : !> \param is_fixed ...
135 : !> \param vcom ...
136 : !> \param ecom ...
137 : !> \par History
138 : !> 2007-11-6: created
139 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
140 : ! **************************************************************************************************
141 3116 : SUBROUTINE compute_vcom(part, is_fixed, vcom, ecom)
142 : TYPE(particle_type), DIMENSION(:), POINTER :: part
143 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
144 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vcom
145 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ecom
146 :
147 : INTEGER :: i
148 : REAL(KIND=dp) :: denom, mass
149 : TYPE(atomic_kind_type), POINTER :: atomic_kind
150 :
151 3116 : vcom = 0.0_dp
152 3116 : denom = 0.0_dp
153 746461 : DO i = 1, SIZE(part)
154 743345 : atomic_kind => part(i)%atomic_kind
155 743345 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
156 746461 : IF (mass /= 0.0) THEN
157 1479322 : SELECT CASE (is_fixed(i))
158 : CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
159 737303 : vcom(1) = vcom(1) + part(i)%v(1)*mass
160 737303 : vcom(2) = vcom(2) + part(i)%v(2)*mass
161 737303 : vcom(3) = vcom(3) + part(i)%v(3)*mass
162 742019 : denom = denom + mass
163 : END SELECT
164 : END IF
165 : END DO
166 12464 : vcom = vcom/denom
167 3116 : IF (PRESENT(ecom)) THEN
168 4364 : ecom = 0.5_dp*denom*SUM(vcom*vcom)
169 : END IF
170 :
171 3116 : END SUBROUTINE compute_vcom
172 :
173 : ! **************************************************************************************************
174 : !> \brief Copy atom velocities into core and shell velocities
175 : !> *** is only used by initialize_velocities below ***
176 : !> \param part ...
177 : !> \param shell_part ...
178 : !> \param core_part ...
179 : !> \par History
180 : !> 2007-11-6: created
181 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
182 : ! **************************************************************************************************
183 8 : SUBROUTINE clone_core_shell_vel(part, shell_part, core_part)
184 : TYPE(particle_type), DIMENSION(:), POINTER :: part, shell_part, core_part
185 :
186 : INTEGER :: i
187 : LOGICAL :: is_shell
188 : TYPE(atomic_kind_type), POINTER :: atomic_kind
189 :
190 776 : DO i = 1, SIZE(part)
191 768 : atomic_kind => part(i)%atomic_kind
192 768 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_shell)
193 776 : IF (is_shell) THEN
194 6144 : shell_part(part(i)%shell_index)%v(:) = part(i)%v(:)
195 6144 : core_part(part(i)%shell_index)%v(:) = part(i)%v(:)
196 : END IF
197 : END DO
198 :
199 8 : END SUBROUTINE clone_core_shell_vel
200 :
201 : ! **************************************************************************************************
202 : !> \brief Compute the kinetic energy. Does not subtract the center of mass kinetic
203 : !> energy.
204 : !> *** is only used by initialize_velocities below ***
205 : !> \param part ...
206 : !> \param ireg ...
207 : !> \return ...
208 : !> \par History
209 : !> 2007-11-6: created
210 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
211 : ! **************************************************************************************************
212 3382 : FUNCTION compute_ekin(part, ireg) RESULT(ekin)
213 : TYPE(particle_type), DIMENSION(:), POINTER :: part
214 : INTEGER, INTENT(IN), OPTIONAL :: ireg
215 : REAL(KIND=dp) :: ekin
216 :
217 : INTEGER :: i
218 : REAL(KIND=dp) :: mass
219 : TYPE(atomic_kind_type), POINTER :: atomic_kind
220 :
221 3382 : NULLIFY (atomic_kind)
222 3382 : ekin = 0.0_dp
223 3382 : IF (PRESENT(ireg)) THEN
224 13756 : DO i = 1, SIZE(part)
225 13756 : IF (part(i)%t_region_index == ireg) THEN
226 4560 : atomic_kind => part(i)%atomic_kind
227 4560 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
228 18240 : ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
229 : END IF
230 : END DO
231 : ELSE
232 746267 : DO i = 1, SIZE(part)
233 743153 : atomic_kind => part(i)%atomic_kind
234 743153 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
235 2975726 : ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
236 : END DO
237 : END IF
238 :
239 3382 : END FUNCTION compute_ekin
240 :
241 : ! **************************************************************************************************
242 : !> \brief Rescale the velocity to mimic the given external kinetic temperature.
243 : !> Optionally also rescale vcom.
244 : !> *** is only used by initialize_velocities below ***
245 : !> \param part ...
246 : !> \param simpar ...
247 : !> \param ekin ...
248 : !> \param vcom ...
249 : !> \param ireg ...
250 : !> \param nfree ...
251 : !> \param temp ...
252 : !> \par History
253 : !> 2007-11-6: created
254 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
255 : ! **************************************************************************************************
256 2053 : SUBROUTINE rescale_vel(part, simpar, ekin, vcom, ireg, nfree, temp)
257 : TYPE(particle_type), DIMENSION(:), POINTER :: part
258 : TYPE(simpar_type), POINTER :: simpar
259 : REAL(KIND=dp), INTENT(INOUT) :: ekin
260 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
261 : OPTIONAL :: vcom
262 : INTEGER, INTENT(IN), OPTIONAL :: ireg, nfree
263 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: temp
264 :
265 : INTEGER :: i, my_ireg, my_nfree
266 : REAL(KIND=dp) :: factor, my_temp
267 :
268 2053 : IF (PRESENT(ireg) .AND. PRESENT(nfree) .AND. PRESENT(temp)) THEN
269 6 : my_ireg = ireg
270 6 : my_nfree = nfree
271 6 : my_temp = temp
272 2047 : ELSEIF (PRESENT(nfree)) THEN
273 0 : my_ireg = 0
274 0 : my_nfree = nfree
275 0 : my_temp = simpar%temp_ext
276 : ELSE
277 2047 : my_ireg = 0
278 2047 : my_nfree = simpar%nfree
279 2047 : my_temp = simpar%temp_ext
280 : END IF
281 2053 : IF (my_nfree /= 0) THEN
282 2041 : factor = my_temp/(2.0_dp*ekin)*REAL(my_nfree, KIND=dp)
283 : ELSE
284 : factor = 0.0_dp
285 : END IF
286 : ! Note:
287 : ! this rescaling is still wrong, it should take the masses into account
288 : ! rescaling is generally not correct, so needs fixing
289 2053 : ekin = ekin*factor
290 2053 : factor = SQRT(factor)
291 2053 : IF (PRESENT(ireg)) THEN
292 582 : DO i = 1, SIZE(part)
293 1158 : IF (part(i)%t_region_index == my_ireg) part(i)%v(:) = factor*part(i)%v(:)
294 : END DO
295 : ELSE
296 444599 : DO i = 1, SIZE(part)
297 1772255 : part(i)%v(:) = factor*part(i)%v(:)
298 : END DO
299 2047 : IF (PRESENT(vcom)) THEN
300 96 : vcom = factor*vcom
301 : END IF
302 : END IF
303 :
304 2053 : END SUBROUTINE rescale_vel
305 :
306 : ! **************************************************************************************************
307 : !> \brief Rescale the velocity of separated regions independently
308 : !> \param part ...
309 : !> \param md_env ...
310 : !> \param simpar ...
311 : !> \par History
312 : !> 2008-11
313 : !> \author MI
314 : ! **************************************************************************************************
315 2 : SUBROUTINE rescale_vel_region(part, md_env, simpar)
316 :
317 : TYPE(particle_type), DIMENSION(:), POINTER :: part
318 : TYPE(md_environment_type), POINTER :: md_env
319 : TYPE(simpar_type), POINTER :: simpar
320 :
321 : INTEGER :: ireg, nfree, nfree0, nfree_done
322 : REAL(KIND=dp) :: ekin, temp
323 : TYPE(thermal_region_type), POINTER :: t_region
324 : TYPE(thermal_regions_type), POINTER :: thermal_regions
325 :
326 2 : NULLIFY (thermal_regions, t_region)
327 :
328 2 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
329 2 : nfree_done = 0
330 6 : DO ireg = 1, thermal_regions%nregions
331 4 : NULLIFY (t_region)
332 4 : t_region => thermal_regions%thermal_region(ireg)
333 4 : nfree = t_region%npart*3
334 4 : ekin = compute_ekin(part, ireg)
335 4 : temp = t_region%temp_expected
336 4 : CALL rescale_vel(part, simpar, ekin, ireg=ireg, nfree=nfree, temp=temp)
337 4 : nfree_done = nfree_done + nfree
338 4 : ekin = compute_ekin(part, ireg)
339 4 : temp = 2.0_dp*ekin/REAL(nfree, dp)*kelvin
340 6 : t_region%temperature = temp
341 : END DO
342 2 : nfree0 = simpar%nfree - nfree_done
343 2 : IF (nfree0 > 0) THEN
344 2 : ekin = compute_ekin(part, 0)
345 2 : CALL rescale_vel(part, simpar, ekin, ireg=0, nfree=nfree0, temp=simpar%temp_ext)
346 2 : ekin = compute_ekin(part, 0)
347 2 : temp = 2.0_dp*ekin/REAL(nfree0, dp)*kelvin
348 2 : thermal_regions%temp_reg0 = temp
349 : END IF
350 2 : END SUBROUTINE rescale_vel_region
351 :
352 : ! **************************************************************************************************
353 : !> \brief subtract center of mass velocity
354 : !> *** is only used by initialize_velocities below ***
355 : !> \param part ...
356 : !> \param is_fixed ...
357 : !> \param vcom ...
358 : !> \par History
359 : !> 2007-11-6: created
360 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
361 : ! **************************************************************************************************
362 2025 : SUBROUTINE subtract_vcom(part, is_fixed, vcom)
363 : TYPE(particle_type), DIMENSION(:), POINTER :: part
364 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
365 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vcom
366 :
367 : INTEGER :: i
368 :
369 442899 : DO i = 1, SIZE(part)
370 2025 : SELECT CASE (is_fixed(i))
371 : CASE (use_perd_x)
372 512 : part(i)%v(2) = part(i)%v(2) - vcom(2)
373 512 : part(i)%v(3) = part(i)%v(3) - vcom(3)
374 : CASE (use_perd_y)
375 512 : part(i)%v(1) = part(i)%v(1) - vcom(1)
376 512 : part(i)%v(3) = part(i)%v(3) - vcom(3)
377 : CASE (use_perd_z)
378 512 : part(i)%v(1) = part(i)%v(1) - vcom(1)
379 512 : part(i)%v(2) = part(i)%v(2) - vcom(2)
380 : CASE (use_perd_xy)
381 512 : part(i)%v(3) = part(i)%v(3) - vcom(3)
382 : CASE (use_perd_xz)
383 0 : part(i)%v(2) = part(i)%v(2) - vcom(2)
384 : CASE (use_perd_yz)
385 0 : part(i)%v(1) = part(i)%v(1) - vcom(1)
386 : CASE (use_perd_none)
387 1747998 : part(i)%v(:) = part(i)%v(:) - vcom(:)
388 : END SELECT
389 : END DO
390 2025 : END SUBROUTINE subtract_vcom
391 :
392 : ! **************************************************************************************************
393 : !> \brief compute the angular velocity
394 : !> *** is only used by initialize_velocities below ***
395 : !> \param part ...
396 : !> \param is_fixed ...
397 : !> \param rcom ...
398 : !> \param vang ...
399 : !> \par History
400 : !> 2007-11-9: created
401 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
402 : ! **************************************************************************************************
403 110 : SUBROUTINE compute_vang(part, is_fixed, rcom, vang)
404 : TYPE(particle_type), DIMENSION(:), POINTER :: part
405 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
406 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rcom
407 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vang
408 :
409 : INTEGER :: i
410 : REAL(KIND=dp) :: mass, proj
411 : REAL(KIND=dp), DIMENSION(3) :: evals, mang, r
412 : REAL(KIND=dp), DIMENSION(3, 3) :: iner
413 : TYPE(atomic_kind_type), POINTER :: atomic_kind
414 :
415 110 : NULLIFY (atomic_kind)
416 110 : mang(:) = 0.0_dp
417 110 : iner(:, :) = 0.0_dp
418 1299 : DO i = 1, SIZE(part)
419 : ! compute angular momentum and inertia tensor
420 110 : SELECT CASE (is_fixed(i))
421 : CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
422 4756 : r(:) = part(i)%r(:) - rcom(:)
423 1189 : atomic_kind => part(i)%atomic_kind
424 1189 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
425 1189 : mang(1) = mang(1) + mass*(r(2)*part(i)%v(3) - r(3)*part(i)%v(2))
426 1189 : mang(2) = mang(2) + mass*(r(3)*part(i)%v(1) - r(1)*part(i)%v(3))
427 1189 : mang(3) = mang(3) + mass*(r(1)*part(i)%v(2) - r(2)*part(i)%v(1))
428 :
429 1189 : iner(1, 1) = iner(1, 1) + mass*(r(2)*r(2) + r(3)*r(3))
430 1189 : iner(2, 2) = iner(2, 2) + mass*(r(3)*r(3) + r(1)*r(1))
431 1189 : iner(3, 3) = iner(3, 3) + mass*(r(1)*r(1) + r(2)*r(2))
432 :
433 1189 : iner(1, 2) = iner(1, 2) - mass*r(1)*r(2)
434 1189 : iner(2, 3) = iner(2, 3) - mass*r(2)*r(3)
435 2378 : iner(3, 1) = iner(3, 1) - mass*r(3)*r(1)
436 : END SELECT
437 : END DO
438 110 : iner(2, 1) = iner(1, 2)
439 110 : iner(3, 2) = iner(2, 3)
440 110 : iner(1, 3) = iner(3, 1)
441 :
442 : ! Take the safest route, i.e. diagonalize the inertia tensor and solve
443 : ! the angular velocity only with the non-zero eigenvalues. A plain inversion
444 : ! would fail for linear molecules.
445 110 : CALL diamat_all(iner, evals)
446 :
447 110 : vang(:) = 0.0_dp
448 440 : DO i = 1, 3
449 440 : IF (evals(i) > 0.0_dp) THEN
450 1276 : proj = SUM(iner(:, i)*mang)/evals(i)
451 319 : vang(1) = vang(1) + proj*iner(1, i)
452 319 : vang(2) = vang(2) + proj*iner(2, i)
453 319 : vang(3) = vang(3) + proj*iner(3, i)
454 : END IF
455 : END DO
456 :
457 110 : END SUBROUTINE compute_vang
458 :
459 : ! **************************************************************************************************
460 : !> \brief subtract the angular velocity
461 : !> *** is only used by initialize_velocities below ***
462 : !> \param part ...
463 : !> \param is_fixed ...
464 : !> \param rcom ...
465 : !> \param vang ...
466 : !> \par History
467 : !> 2007-11-9: created
468 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
469 : ! **************************************************************************************************
470 6 : SUBROUTINE subtract_vang(part, is_fixed, rcom, vang)
471 : TYPE(particle_type), DIMENSION(:), POINTER :: part
472 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
473 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rcom, vang
474 :
475 : INTEGER :: i
476 : REAL(KIND=dp), DIMENSION(3) :: r
477 :
478 52 : DO i = 1, SIZE(part)
479 184 : r(:) = part(i)%r(:) - rcom(:)
480 6 : SELECT CASE (is_fixed(i))
481 : CASE (use_perd_x)
482 0 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
483 0 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
484 : CASE (use_perd_y)
485 0 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
486 0 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
487 : CASE (use_perd_z)
488 0 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
489 0 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
490 : CASE (use_perd_xy)
491 0 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
492 : CASE (use_perd_xz)
493 0 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
494 : CASE (use_perd_yz)
495 0 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
496 : CASE (use_perd_none)
497 46 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
498 46 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
499 46 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
500 : END SELECT
501 : END DO
502 :
503 6 : END SUBROUTINE subtract_vang
504 :
505 : ! **************************************************************************************************
506 : !> \brief Initializes the velocities to the Maxwellian distribution
507 : !> \param simpar ...
508 : !> \param part ...
509 : !> \param force_env ...
510 : !> \param globenv ...
511 : !> \param md_env ...
512 : !> \param molecule_kinds ...
513 : !> \param label ...
514 : !> \param print_section ...
515 : !> \param subsys_section ...
516 : !> \param shell_present ...
517 : !> \param shell_part ...
518 : !> \param core_part ...
519 : !> \param force_rescaling ...
520 : !> \param para_env ...
521 : !> \param write_binary_restart_file ...
522 : !> \par History
523 : !> - is_fixed removed from particle_type
524 : !> - 2007-11-07: Cleanup (TV)
525 : !> - 2007-11-09: Added angvel_zero feature
526 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>
527 : ! **************************************************************************************************
528 1791 : SUBROUTINE initialize_velocities(simpar, &
529 : part, &
530 : force_env, &
531 : globenv, &
532 : md_env, &
533 : molecule_kinds, &
534 : label, &
535 : print_section, &
536 : subsys_section, &
537 : shell_present, &
538 : shell_part, &
539 : core_part, &
540 : force_rescaling, &
541 : para_env, &
542 : write_binary_restart_file)
543 :
544 : TYPE(simpar_type), POINTER :: simpar
545 : TYPE(particle_type), DIMENSION(:), POINTER :: part
546 : TYPE(force_env_type), POINTER :: force_env
547 : TYPE(global_environment_type), POINTER :: globenv
548 : TYPE(md_environment_type), POINTER :: md_env
549 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
550 : CHARACTER(LEN=*), INTENT(IN) :: label
551 : TYPE(section_vals_type), POINTER :: print_section, subsys_section
552 : LOGICAL, INTENT(IN) :: shell_present
553 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
554 : LOGICAL, INTENT(IN) :: force_rescaling
555 : TYPE(mp_para_env_type), POINTER :: para_env
556 : LOGICAL, INTENT(IN) :: write_binary_restart_file
557 :
558 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_velocities'
559 :
560 : INTEGER :: handle, i, ifixd, imolecule_kind, iw, &
561 : natoms
562 : INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
563 : LOGICAL :: success
564 : REAL(KIND=dp) :: ecom, ekin, mass, mass_tot, temp, tmp_r1
565 : REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
566 : TYPE(atomic_kind_type), POINTER :: atomic_kind
567 : TYPE(cell_type), POINTER :: cell
568 : TYPE(cp_logger_type), POINTER :: logger
569 1791 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
570 1791 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
571 : TYPE(molecule_kind_type), POINTER :: molecule_kind
572 : TYPE(section_vals_type), POINTER :: md_section, root_section, vib_section
573 :
574 1791 : CALL timeset(routineN, handle)
575 :
576 : ! Initializing parameters
577 1791 : natoms = SIZE(part)
578 1791 : NULLIFY (atomic_kind, fixd_list, logger, molecule_kind)
579 1791 : NULLIFY (molecule_kind_set)
580 :
581 : ! Logging
582 1791 : logger => cp_get_default_logger()
583 1791 : iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
584 :
585 : ! Build a list of all fixed atoms (if any)
586 5373 : ALLOCATE (is_fixed(natoms))
587 :
588 1791 : is_fixed = use_perd_none
589 1791 : molecule_kind_set => molecule_kinds%els
590 66581 : DO imolecule_kind = 1, molecule_kinds%n_els
591 64790 : molecule_kind => molecule_kind_set(imolecule_kind)
592 64790 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
593 66581 : IF (ASSOCIATED(fixd_list)) THEN
594 5432 : DO ifixd = 1, SIZE(fixd_list)
595 5432 : IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
596 : END DO
597 : END IF
598 : END DO
599 :
600 : ! Compute the total mass when needed
601 1791 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
602 : simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
603 : mass_tot = 0.0_dp
604 1006 : DO i = 1, natoms
605 1000 : atomic_kind => part(i)%atomic_kind
606 1000 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
607 1006 : mass_tot = mass_tot + mass
608 : END DO
609 6 : simpar%v_shock = simpar%v_shock*SQRT(mass_tot)
610 : END IF
611 :
612 : CALL read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
613 1791 : shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
614 1791 : IF (.NOT. success) THEN
615 3052 : SELECT CASE (simpar%initialization_method)
616 : CASE (md_init_default)
617 : CALL generate_velocities(simpar, part, force_env, globenv, md_env, shell_present, &
618 1525 : shell_part, core_part, is_fixed, iw)
619 : CASE (md_init_vib)
620 2 : CALL force_env_get(force_env=force_env, root_section=root_section)
621 2 : md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
622 2 : vib_section => section_vals_get_subs_vals(root_section, "VIBRATIONAL_ANALYSIS")
623 : CALL generate_coords_vels_vib(simpar, &
624 : part, &
625 : md_section, &
626 : vib_section, &
627 : force_env, &
628 : globenv, &
629 : shell_present, &
630 : shell_part, &
631 : core_part, &
632 2 : is_fixed)
633 : ! update restart file for the modified coordinates and velocities
634 1529 : CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
635 : END SELECT
636 : END IF
637 :
638 1791 : IF (iw > 0) THEN
639 : WRITE (iw, '(/,T2,A)') &
640 826 : 'MD_VEL| '//TRIM(ADJUSTL(label))
641 : ! Recompute vcom, ecom and ekin for IO
642 826 : CALL compute_vcom(part, is_fixed, vcom, ecom)
643 826 : ekin = compute_ekin(part) - ecom
644 826 : IF (simpar%nfree == 0) THEN
645 6 : CPASSERT(ekin == 0.0_dp)
646 6 : temp = 0.0_dp
647 : ELSE
648 820 : temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
649 : END IF
650 826 : tmp_r1 = cp_unit_from_cp2k(temp, "K")
651 : WRITE (iw, '(T2,A,T61,F20.6)') &
652 826 : 'MD_VEL| Initial temperature [K]', tmp_r1
653 : WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
654 826 : 'MD_VEL| COM velocity', vcom(1:3)
655 :
656 : ! compute and log rcom and vang if not periodic
657 826 : CALL force_env_get(force_env, cell=cell)
658 3304 : IF (SUM(cell%perd(1:3)) == 0) THEN
659 64 : CALL compute_rcom(part, is_fixed, rcom)
660 64 : CALL compute_vang(part, is_fixed, rcom, vang)
661 : WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
662 64 : 'MD_VEL| COM position', rcom(1:3)
663 : WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
664 64 : 'MD_VEL| Angular velocity', vang(1:3)
665 : END IF
666 : END IF
667 :
668 1791 : DEALLOCATE (is_fixed)
669 1791 : CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
670 1791 : CALL timestop(handle)
671 :
672 3582 : END SUBROUTINE initialize_velocities
673 :
674 : ! **************************************************************************************************
675 : !> \brief Read velocities from binary restart file if available
676 : !> \param simpar ...
677 : !> \param part ...
678 : !> \param force_env ...
679 : !> \param md_env ...
680 : !> \param subsys_section ...
681 : !> \param shell_present ...
682 : !> \param shell_part ...
683 : !> \param core_part ...
684 : !> \param force_rescaling ...
685 : !> \param para_env ...
686 : !> \param is_fixed ...
687 : !> \param success ...
688 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>
689 : ! **************************************************************************************************
690 12537 : SUBROUTINE read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
691 1791 : shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
692 : TYPE(simpar_type), POINTER :: simpar
693 : TYPE(particle_type), DIMENSION(:), POINTER :: part
694 : TYPE(force_env_type), POINTER :: force_env
695 : TYPE(md_environment_type), POINTER :: md_env
696 : TYPE(section_vals_type), POINTER :: subsys_section
697 : LOGICAL, INTENT(IN) :: shell_present
698 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
699 : LOGICAL, INTENT(IN) :: force_rescaling
700 : TYPE(mp_para_env_type), POINTER :: para_env
701 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
702 : LOGICAL, INTENT(OUT) :: success
703 :
704 : INTEGER :: i, natoms, nshell, shell_index
705 : LOGICAL :: atomvel_explicit, atomvel_read, corevel_explicit, corevel_read, is_ok, &
706 : rescale_regions, shellvel_explicit, shellvel_read
707 : REAL(KIND=dp) :: ecom, ekin, fac_massc, fac_masss, mass
708 : REAL(KIND=dp), DIMENSION(3) :: vc, vcom, vs
709 1791 : REAL(KIND=dp), DIMENSION(:), POINTER :: vel
710 : TYPE(atomic_kind_type), POINTER :: atomic_kind
711 : TYPE(cell_type), POINTER :: cell
712 : TYPE(cp_sll_val_type), POINTER :: atom_list, core_list, shell_list
713 : TYPE(section_vals_type), POINTER :: atomvel_section, corevel_section, &
714 : shellvel_section
715 : TYPE(shell_kind_type), POINTER :: shell
716 : TYPE(thermal_regions_type), POINTER :: thermal_regions
717 : TYPE(val_type), POINTER :: val
718 :
719 : ! Initializing parameters
720 :
721 1791 : success = .FALSE.
722 1791 : natoms = SIZE(part)
723 : atomvel_read = .FALSE.
724 : corevel_read = .FALSE.
725 : shellvel_read = .FALSE.
726 1791 : NULLIFY (vel, atomic_kind, atom_list, core_list, shell_list)
727 1791 : NULLIFY (atomvel_section, shellvel_section, corevel_section)
728 1791 : NULLIFY (cell, shell, thermal_regions, val)
729 1791 : CALL force_env_get(force_env, cell=cell)
730 :
731 : ! Core-Shell Model
732 1791 : nshell = 0
733 1791 : IF (shell_present) THEN
734 132 : CPASSERT(ASSOCIATED(core_part))
735 132 : CPASSERT(ASSOCIATED(shell_part))
736 132 : nshell = SIZE(shell_part)
737 : END IF
738 :
739 1791 : atomvel_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
740 1791 : shellvel_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
741 1791 : corevel_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
742 :
743 : ! Read or initialize the particle velocities
744 1791 : CALL section_vals_get(atomvel_section, explicit=atomvel_explicit)
745 1791 : CALL section_vals_get(shellvel_section, explicit=shellvel_explicit)
746 1791 : CALL section_vals_get(corevel_section, explicit=corevel_explicit)
747 1791 : CPASSERT(shellvel_explicit .EQV. corevel_explicit)
748 :
749 : CALL read_binary_velocities("", part, force_env%root_section, para_env, &
750 1791 : subsys_section, atomvel_read, cell)
751 : CALL read_binary_velocities("SHELL", shell_part, force_env%root_section, para_env, &
752 1791 : subsys_section, shellvel_read, cell)
753 : CALL read_binary_velocities("CORE", core_part, force_env%root_section, para_env, &
754 1791 : subsys_section, corevel_read, cell)
755 :
756 1791 : IF (.NOT. (atomvel_explicit .OR. atomvel_read)) RETURN
757 264 : success = .TRUE.
758 :
759 264 : IF (.NOT. atomvel_read) THEN
760 : ! Read the atom velocities if explicitly given in the input file
761 218 : CALL section_vals_list_get(atomvel_section, "_DEFAULT_KEYWORD_", list=atom_list)
762 52296 : DO i = 1, natoms
763 52078 : is_ok = cp_sll_val_next(atom_list, val)
764 52078 : CALL val_get(val, r_vals=vel)
765 364546 : part(i)%v = vel
766 52296 : CALL cell_transform_input_cartesian(cell, part(i)%v)
767 : END DO
768 : END IF
769 56758 : DO i = 1, natoms
770 264 : SELECT CASE (is_fixed(i))
771 : CASE (use_perd_x)
772 0 : part(i)%v(1) = 0.0_dp
773 : CASE (use_perd_y)
774 0 : part(i)%v(2) = 0.0_dp
775 : CASE (use_perd_z)
776 0 : part(i)%v(3) = 0.0_dp
777 : CASE (use_perd_xy)
778 0 : part(i)%v(1) = 0.0_dp
779 0 : part(i)%v(2) = 0.0_dp
780 : CASE (use_perd_xz)
781 0 : part(i)%v(1) = 0.0_dp
782 0 : part(i)%v(3) = 0.0_dp
783 : CASE (use_perd_yz)
784 0 : part(i)%v(2) = 0.0_dp
785 0 : part(i)%v(3) = 0.0_dp
786 : CASE (use_perd_xyz)
787 56572 : part(i)%v = 0.0_dp
788 : END SELECT
789 : END DO
790 264 : IF (shell_present) THEN
791 48 : IF (shellvel_explicit) THEN
792 : ! If the atoms positions are given (?) and core and shell velocities are
793 : ! present in the input, read the latter.
794 24 : CALL section_vals_list_get(shellvel_section, "_DEFAULT_KEYWORD_", list=shell_list)
795 24 : CALL section_vals_list_get(corevel_section, "_DEFAULT_KEYWORD_", list=core_list)
796 2328 : DO i = 1, nshell
797 2304 : is_ok = cp_sll_val_next(shell_list, val)
798 2304 : CALL val_get(val, r_vals=vel)
799 16128 : shell_part(i)%v = vel
800 2304 : CALL cell_transform_input_cartesian(cell, shell_part(i)%v)
801 2304 : is_ok = cp_sll_val_next(core_list, val)
802 2304 : CALL val_get(val, r_vals=vel)
803 16128 : core_part(i)%v = vel
804 2328 : CALL cell_transform_input_cartesian(cell, core_part(i)%v)
805 : END DO
806 : ELSE
807 24 : IF (.NOT. (shellvel_read .AND. corevel_read)) THEN
808 : ! Otherwise, just copy atom velocties into shell and core velocities.
809 8 : CALL clone_core_shell_vel(part, shell_part, core_part)
810 : END IF
811 : END IF
812 : END IF
813 :
814 : ! compute vcom, ecom and ekin
815 264 : CALL compute_vcom(part, is_fixed, vcom, ecom)
816 264 : ekin = compute_ekin(part) - ecom
817 :
818 264 : IF (simpar%do_thermal_region) THEN
819 12 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
820 12 : IF (ASSOCIATED(thermal_regions)) THEN
821 12 : rescale_regions = thermal_regions%force_rescaling
822 : END IF
823 : ELSE
824 : rescale_regions = .FALSE.
825 : END IF
826 264 : IF (simpar%nfree /= 0 .AND. (force_rescaling .OR. rescale_regions)) THEN
827 24 : IF (simpar%do_thermal_region) THEN
828 0 : CALL rescale_vel_region(part, md_env, simpar)
829 : ELSE
830 24 : CALL rescale_vel(part, simpar, ekin, vcom=vcom)
831 : END IF
832 :
833 : ! After rescaling, the core and shell velocities must also adapt.
834 1894 : DO i = 1, natoms
835 1870 : shell_index = part(i)%shell_index
836 2134 : IF (shell_present .AND. shell_index /= 0) THEN
837 0 : atomic_kind => part(i)%atomic_kind
838 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
839 0 : fac_masss = shell%mass_shell/mass
840 0 : fac_massc = shell%mass_core/mass
841 0 : vs = shell_part(shell_index)%v
842 0 : vc = core_part(shell_index)%v
843 :
844 0 : shell_part(shell_index)%v(1) = part(i)%v(1) + fac_massc*(vs(1) - vc(1))
845 0 : shell_part(shell_index)%v(2) = part(i)%v(2) + fac_massc*(vs(2) - vc(2))
846 0 : shell_part(shell_index)%v(3) = part(i)%v(3) + fac_massc*(vs(3) - vc(3))
847 0 : core_part(shell_index)%v(1) = part(i)%v(1) + fac_masss*(vc(1) - vs(1))
848 0 : core_part(shell_index)%v(2) = part(i)%v(2) + fac_masss*(vc(2) - vs(2))
849 0 : core_part(shell_index)%v(3) = part(i)%v(3) + fac_masss*(vc(3) - vs(3))
850 : END IF
851 : END DO
852 : END IF
853 1791 : END SUBROUTINE read_input_velocities
854 :
855 : ! **************************************************************************************************
856 : !> \brief Initializing velocities AND positions randomly on all processors, based on vibrational
857 : !> modes of the system, so that the starting coordinates are already very close to
858 : !> canonical ensumble corresponding to temperature of a head bath.
859 : !> \param simpar : MD simulation parameters
860 : !> \param particles : global array of all particles
861 : !> \param md_section : MD input subsection
862 : !> \param vib_section : vibrational analysis input section
863 : !> \param force_env : force environment container
864 : !> \param global_env : global environment container
865 : !> \param shell_present : if core-shell model is used
866 : !> \param shell_particles : global array of all shell particles in shell model
867 : !> \param core_particles : global array of all core particles in shell model
868 : !> \param is_fixed : array of size of total number of atoms, that determines if any
869 : !> cartesian components are fixed
870 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
871 : ! **************************************************************************************************
872 2 : SUBROUTINE generate_coords_vels_vib(simpar, &
873 : particles, &
874 : md_section, &
875 : vib_section, &
876 : force_env, &
877 : global_env, &
878 : shell_present, &
879 : shell_particles, &
880 : core_particles, &
881 2 : is_fixed)
882 : TYPE(simpar_type), POINTER :: simpar
883 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
884 : TYPE(section_vals_type), POINTER :: md_section, vib_section
885 : TYPE(force_env_type), POINTER :: force_env
886 : TYPE(global_environment_type), POINTER :: global_env
887 : LOGICAL, INTENT(IN) :: shell_present
888 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_particles, core_particles
889 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
890 :
891 : INTEGER :: dof, fixed_dof, iatom, ii, imode, &
892 : my_dof, natoms, shell_index
893 : REAL(KIND=dp) :: Erand, mass, my_phase, ratio, temperature
894 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, phase, random
895 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dr, eigenvectors
896 : TYPE(atomic_kind_type), POINTER :: atomic_kind
897 : TYPE(mp_para_env_type), POINTER :: para_env
898 2 : TYPE(rng_stream_type), ALLOCATABLE :: random_stream
899 :
900 2 : CALL cite_reference(West2006)
901 2 : natoms = SIZE(particles)
902 2 : temperature = simpar%temp_ext
903 2 : my_dof = 3*natoms
904 6 : ALLOCATE (eigenvalues(my_dof))
905 8 : ALLOCATE (eigenvectors(my_dof, my_dof))
906 4 : ALLOCATE (phase(my_dof))
907 4 : ALLOCATE (random(my_dof))
908 6 : ALLOCATE (dr(3, natoms))
909 2 : CALL force_env_get(force_env=force_env, para_env=para_env)
910 : ! read vibration modes
911 : CALL read_vib_eigs_unformatted(md_section, &
912 : vib_section, &
913 : para_env, &
914 : dof, &
915 : eigenvalues, &
916 2 : eigenvectors)
917 2 : IF (my_dof /= dof) THEN
918 : CALL cp_abort(__LOCATION__, &
919 : "number of degrees of freedom in vibrational analysis data "// &
920 0 : "do not match total number of cartesian degrees of freedom")
921 : END IF
922 : ! read phases
923 2 : CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%PHASE", r_val=my_phase)
924 2 : my_phase = MIN(1.0_dp, my_phase)
925 : ! generate random numbers
926 2 : random_stream = rng_stream_type(name="MD_INIT_VIB", distribution_type=UNIFORM)
927 2 : CALL random_stream%fill(random)
928 2 : IF (my_phase < 0.0_dp) THEN
929 0 : CALL random_stream%fill(phase)
930 : ELSE
931 20 : phase = my_phase
932 : END IF
933 2 : DEALLOCATE (random_stream)
934 :
935 : ! the first three modes are acoustic with zero frequencies,
936 : ! exclude these from considerations
937 2 : my_dof = dof - 3
938 : ! randomly selects energy from distribution about kT, all
939 : ! energies are scaled so that the sum over vibration modes gives
940 : ! exactly my_dof*kT. Note that k = 1.0 in atomic units
941 2 : Erand = 0.0_dp
942 14 : DO imode = 4, dof
943 14 : Erand = Erand - temperature*LOG(1.0_dp - random(imode))
944 : END DO
945 : ! need to take into account of fixed constraints too
946 2 : fixed_dof = 0
947 8 : DO iatom = 1, natoms
948 2 : SELECT CASE (is_fixed(iatom))
949 : CASE (use_perd_x, use_perd_y, use_perd_z)
950 0 : fixed_dof = fixed_dof + 1
951 : CASE (use_perd_xy, use_perd_xz, use_perd_yz)
952 0 : fixed_dof = fixed_dof + 2
953 : CASE (use_perd_xyz)
954 6 : fixed_dof = fixed_dof + 3
955 : END SELECT
956 : END DO
957 2 : my_dof = my_dof - fixed_dof
958 2 : ratio = REAL(my_dof, KIND=dp)*temperature/Erand
959 : ! update velocities AND positions
960 8 : DO iatom = 1, natoms
961 6 : atomic_kind => particles(iatom)%atomic_kind
962 6 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
963 8 : SELECT CASE (is_fixed(iatom))
964 : CASE (use_perd_x)
965 0 : DO ii = 2, 3
966 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
967 0 : eigenvectors, random, phase, dof, ratio)
968 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
969 : eigenvectors, random, phase, dof, &
970 0 : ratio)
971 : END DO
972 : CASE (use_perd_y)
973 0 : DO ii = 1, 3, 2
974 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
975 0 : eigenvectors, random, phase, dof, ratio)
976 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
977 : eigenvectors, random, phase, dof, &
978 0 : ratio)
979 : END DO
980 : CASE (use_perd_z)
981 0 : DO ii = 1, 2
982 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
983 0 : eigenvectors, random, phase, dof, ratio)
984 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
985 : eigenvectors, random, phase, dof, &
986 0 : ratio)
987 : END DO
988 : CASE (use_perd_xy)
989 : dr(3, iatom) = dr_from_vib_data(iatom, 3, mass, temperature, eigenvalues, &
990 0 : eigenvectors, random, phase, dof, ratio)
991 : particles(iatom)%v(3) = dv_from_vib_data(iatom, 3, mass, temperature, &
992 : eigenvectors, random, phase, dof, &
993 0 : ratio)
994 : CASE (use_perd_xz)
995 : dr(2, iatom) = dr_from_vib_data(iatom, 2, mass, temperature, eigenvalues, &
996 0 : eigenvectors, random, phase, dof, ratio)
997 : particles(iatom)%v(2) = dv_from_vib_data(iatom, 2, mass, temperature, &
998 : eigenvectors, random, phase, dof, &
999 0 : ratio)
1000 : CASE (use_perd_yz)
1001 : dr(1, iatom) = dr_from_vib_data(iatom, 1, mass, temperature, eigenvalues, &
1002 0 : eigenvectors, random, phase, dof, ratio)
1003 : particles(iatom)%v(1) = dv_from_vib_data(iatom, 1, mass, temperature, &
1004 : eigenvectors, random, phase, dof, &
1005 0 : ratio)
1006 : CASE (use_perd_none)
1007 24 : DO ii = 1, 3
1008 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
1009 18 : eigenvectors, random, phase, dof, ratio)
1010 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
1011 : eigenvectors, random, phase, dof, &
1012 24 : ratio)
1013 : END DO
1014 : END SELECT
1015 : END DO ! iatom
1016 : ! free memory
1017 2 : DEALLOCATE (eigenvalues)
1018 2 : DEALLOCATE (eigenvectors)
1019 2 : DEALLOCATE (phase)
1020 2 : DEALLOCATE (random)
1021 : ! update particle coordinates
1022 8 : DO iatom = 1, natoms
1023 26 : particles(iatom)%r(:) = particles(iatom)%r(:) + dr(:, iatom)
1024 : END DO
1025 : ! update core-shell model coordinates
1026 2 : IF (shell_present) THEN
1027 : ! particles have moved, and for core-shell model this means
1028 : ! the cores and shells must also move by the same amount. The
1029 : ! shell positions will then be optimised if needed
1030 0 : shell_index = particles(iatom)%shell_index
1031 0 : IF (shell_index /= 0) THEN
1032 : core_particles(shell_index)%r(:) = core_particles(shell_index)%r(:) + &
1033 0 : dr(:, iatom)
1034 : shell_particles(shell_index)%r(:) = shell_particles(shell_index)%r(:) + &
1035 0 : dr(:, iatom)
1036 : END IF
1037 : CALL optimize_shell_core(force_env, &
1038 : particles, &
1039 : shell_particles, &
1040 : core_particles, &
1041 0 : global_env)
1042 : END IF
1043 : ! cleanup
1044 2 : DEALLOCATE (dr)
1045 4 : END SUBROUTINE generate_coords_vels_vib
1046 :
1047 : ! **************************************************************************************************
1048 : !> \brief calculates componbent of initial velocity of an atom from vibreational modes
1049 : !> \param iatom : global atomic index
1050 : !> \param icart : cartesian index (1, 2 or 3)
1051 : !> \param mass : atomic mass
1052 : !> \param temperature : target temperature of canonical ensemble
1053 : !> \param eigenvalues : array containing all cartesian vibrational mode eigenvalues (frequencies)
1054 : !> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
1055 : !> (displacements)
1056 : !> \param random : array containing uniform distributed random numbers, must be the size
1057 : !> of 3*natoms. Numbers must be between 0 and 1
1058 : !> \param phase : array containing numbers between 0 and 1 that determines for each
1059 : !> vibration mode the ratio of potential energy vs kinetic energy contribution
1060 : !> to total energy
1061 : !> \param dof : total number of degrees of freedom, = 3*natoms
1062 : !> \param scale : scale to make sure the sum of vibrational modes give the correct energy
1063 : !> \return : outputs icart-th cartesian component of initial position of atom iatom
1064 : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
1065 : ! **************************************************************************************************
1066 18 : PURE FUNCTION dr_from_vib_data(iatom, &
1067 : icart, &
1068 : mass, &
1069 : temperature, &
1070 36 : eigenvalues, &
1071 18 : eigenvectors, &
1072 18 : random, &
1073 18 : phase, &
1074 : dof, &
1075 : scale) &
1076 : RESULT(res)
1077 : INTEGER, INTENT(IN) :: iatom, icart
1078 : REAL(KIND=dp), INTENT(IN) :: mass, temperature
1079 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
1080 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvectors
1081 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: random, phase
1082 : INTEGER, INTENT(IN) :: dof
1083 : REAL(KIND=dp), INTENT(IN) :: scale
1084 : REAL(KIND=dp) :: res
1085 :
1086 : INTEGER :: imode, ind
1087 :
1088 18 : res = 0.0_dp
1089 : ! assuming the eigenvalues are sorted in ascending order, the
1090 : ! first three modes are acoustic with zero frequencies. These are
1091 : ! excluded from considerations, and should have been reflected in
1092 : ! the calculation of scale outside this function
1093 18 : IF (mass > 0.0_dp) THEN
1094 : ! eigenvector rows assumed to be grouped in atomic blocks
1095 18 : ind = (iatom - 1)*3 + icart
1096 126 : DO imode = 4, dof
1097 : res = res + &
1098 : SQRT(-2.0_dp*scale*temperature*LOG(1 - random(imode))/mass)/ &
1099 : eigenvalues(imode)* &
1100 : eigenvectors(ind, imode)* &
1101 126 : COS(2.0_dp*pi*phase(imode))
1102 : END DO
1103 : END IF
1104 18 : END FUNCTION dr_from_vib_data
1105 :
1106 : ! **************************************************************************************************
1107 : !> \brief calculates componbent of initial velocity of an atom from vibreational modes
1108 : !> \param iatom : global atomic index
1109 : !> \param icart : cartesian index (1, 2 or 3)
1110 : !> \param mass : atomic mass
1111 : !> \param temperature : target temperature of canonical ensemble
1112 : !> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
1113 : !> (displacements)
1114 : !> \param random : array containing uniform distributed random numbers, must be the size
1115 : !> of 3*natoms. Numbers must be between 0 and 1
1116 : !> \param phase : array containing numbers between 0 and 1 that determines for each
1117 : !> vibration mode the ratio of potential energy vs kinetic energy contribution
1118 : !> to total energy
1119 : !> \param dof : total number of degrees of freedom, = 3*natoms
1120 : !> \param scale : scale to make sure the sum of vibrational modes give the correct energy
1121 : !> \return : outputs icart-th cartesian component of initial velocity of atom iatom
1122 : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
1123 : ! **************************************************************************************************
1124 18 : PURE FUNCTION dv_from_vib_data(iatom, &
1125 : icart, &
1126 : mass, &
1127 : temperature, &
1128 36 : eigenvectors, &
1129 18 : random, &
1130 18 : phase, &
1131 : dof, &
1132 : scale) &
1133 : RESULT(res)
1134 : INTEGER, INTENT(IN) :: iatom, icart
1135 : REAL(KIND=dp), INTENT(IN) :: mass, temperature
1136 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvectors
1137 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: random, phase
1138 : INTEGER, INTENT(IN) :: dof
1139 : REAL(KIND=dp), INTENT(IN) :: scale
1140 : REAL(KIND=dp) :: res
1141 :
1142 : INTEGER :: imode, ind
1143 :
1144 18 : res = 0.0_dp
1145 : ! assuming the eigenvalues are sorted in ascending order, the
1146 : ! first three modes are acoustic with zero frequencies. These are
1147 : ! excluded from considerations, and should have been reflected in
1148 : ! the calculation of scale outside this function
1149 18 : IF (mass > 0.0_dp) THEN
1150 : ! eigenvector rows assumed to be grouped in atomic blocks
1151 18 : ind = (iatom - 1)*3 + icart
1152 126 : DO imode = 4, dof
1153 : res = res - &
1154 : SQRT(-2.0_dp*scale*temperature*LOG(1 - random(imode))/mass)* &
1155 : eigenvectors(ind, imode)* &
1156 126 : SIN(2.0_dp*pi*phase(imode))
1157 : END DO
1158 : END IF
1159 18 : END FUNCTION dv_from_vib_data
1160 :
1161 : ! **************************************************************************************************
1162 : !> \brief Initializing velocities deterministically on all processors, if not given in input
1163 : !> \param simpar ...
1164 : !> \param part ...
1165 : !> \param force_env ...
1166 : !> \param globenv ...
1167 : !> \param md_env ...
1168 : !> \param shell_present ...
1169 : !> \param shell_part ...
1170 : !> \param core_part ...
1171 : !> \param is_fixed ...
1172 : !> \param iw ...
1173 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1174 : ! **************************************************************************************************
1175 1525 : SUBROUTINE generate_velocities(simpar, part, force_env, globenv, md_env, &
1176 1525 : shell_present, shell_part, core_part, is_fixed, iw)
1177 : TYPE(simpar_type), POINTER :: simpar
1178 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1179 : TYPE(force_env_type), POINTER :: force_env
1180 : TYPE(global_environment_type), POINTER :: globenv
1181 : TYPE(md_environment_type), POINTER :: md_env
1182 : LOGICAL, INTENT(IN) :: shell_present
1183 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
1184 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1185 : INTEGER, INTENT(IN) :: iw
1186 :
1187 : INTEGER :: i, natoms
1188 : REAL(KIND=dp) :: mass
1189 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1190 :
1191 1525 : NULLIFY (atomic_kind)
1192 1525 : natoms = SIZE(part)
1193 :
1194 437399 : DO i = 1, natoms
1195 435874 : atomic_kind => part(i)%atomic_kind
1196 435874 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1197 435874 : part(i)%v(1) = 0.0_dp
1198 435874 : part(i)%v(2) = 0.0_dp
1199 435874 : part(i)%v(3) = 0.0_dp
1200 437399 : IF (mass /= 0.0) THEN
1201 435940 : SELECT CASE (is_fixed(i))
1202 : CASE (use_perd_x)
1203 512 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1204 512 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1205 : CASE (use_perd_y)
1206 512 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1207 512 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1208 : CASE (use_perd_z)
1209 512 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1210 512 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1211 : CASE (use_perd_xy)
1212 512 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1213 : CASE (use_perd_xz)
1214 0 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1215 : CASE (use_perd_yz)
1216 0 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1217 : CASE (use_perd_none)
1218 430262 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1219 430262 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1220 865690 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1221 : END SELECT
1222 : END IF
1223 : END DO
1224 :
1225 1525 : CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1226 1525 : CALL soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
1227 :
1228 : ! Initialize the core and the shell velocity. Atom velocities are just
1229 : ! copied so that the initial relative core-shell velocity is zero.
1230 1525 : IF (shell_present) THEN
1231 84 : CALL optimize_shell_core(force_env, part, shell_part, core_part, globenv)
1232 : END IF
1233 1525 : END SUBROUTINE generate_velocities
1234 :
1235 : ! **************************************************************************************************
1236 : !> \brief Direct velocities along a low-curvature direction in order to
1237 : !> favors MD trajectories to cross rapidly over small energy barriers
1238 : !> into neighboring basins.
1239 : !> \param simpar ...
1240 : !> \param part ...
1241 : !> \param force_env ...
1242 : !> \param md_env ...
1243 : !> \param is_fixed ...
1244 : !> \param iw ...
1245 : !> \author Ole Schuett
1246 : ! **************************************************************************************************
1247 1525 : SUBROUTINE soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
1248 : TYPE(simpar_type), POINTER :: simpar
1249 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1250 : TYPE(force_env_type), POINTER :: force_env
1251 : TYPE(md_environment_type), POINTER :: md_env
1252 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1253 : INTEGER, INTENT(IN) :: iw
1254 :
1255 : INTEGER :: i, k
1256 1500 : REAL(KIND=dp), DIMENSION(SIZE(part), 3) :: F, F_t, N, x0
1257 :
1258 1525 : IF (simpar%soften_nsteps <= 0) RETURN !nothing todo
1259 :
1260 275 : IF (ANY(is_fixed /= use_perd_none)) &
1261 0 : CPABORT("Velocitiy softening with constraints is not supported.")
1262 :
1263 : !backup positions
1264 275 : DO i = 1, SIZE(part)
1265 1025 : x0(i, :) = part(i)%r
1266 : END DO
1267 :
1268 525 : DO k = 1, simpar%soften_nsteps
1269 :
1270 : !use normalized velocities as displace direction
1271 5500 : DO i = 1, SIZE(part)
1272 20500 : N(i, :) = part(i)%v
1273 : END DO
1274 33500 : N = N/SQRT(SUM(N**2))
1275 :
1276 : ! displace system temporarly to calculate forces
1277 5500 : DO i = 1, SIZE(part)
1278 20500 : part(i)%r = part(i)%r + simpar%soften_delta*N(i, :)
1279 : END DO
1280 500 : CALL force_env_calc_energy_force(force_env)
1281 :
1282 : ! calculate velocity update direction F_t
1283 5500 : DO i = 1, SIZE(part)
1284 20500 : F(i, :) = part(i)%f
1285 : END DO
1286 33500 : F_t = F - N*SUM(N*F)
1287 :
1288 : ! restore positions and update velocities
1289 5500 : DO i = 1, SIZE(part)
1290 20000 : part(i)%r = x0(i, :)
1291 20500 : part(i)%v = part(i)%v + simpar%soften_alpha*F_t(i, :)
1292 : END DO
1293 :
1294 525 : CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1295 : END DO
1296 :
1297 25 : IF (iw > 0) THEN
1298 0 : WRITE (iw, "(A,T71, I10)") " Velocities softening Steps: ", simpar%soften_nsteps
1299 0 : WRITE (iw, "(A,T71, E10.3)") " Velocities softening NORM(F_t): ", SQRT(SUM(F_t**2))
1300 : END IF
1301 25 : END SUBROUTINE soften_velocities
1302 :
1303 : ! **************************************************************************************************
1304 : !> \brief Scale velocities according to temperature and remove rigid body motion.
1305 : !> \param simpar ...
1306 : !> \param part ...
1307 : !> \param force_env ...
1308 : !> \param md_env ...
1309 : !> \param is_fixed ...
1310 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1311 : ! **************************************************************************************************
1312 2025 : SUBROUTINE normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1313 : TYPE(simpar_type), POINTER :: simpar
1314 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1315 : TYPE(force_env_type), POINTER :: force_env
1316 : TYPE(md_environment_type), POINTER :: md_env
1317 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1318 :
1319 : REAL(KIND=dp) :: ekin
1320 : REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
1321 : TYPE(cell_type), POINTER :: cell
1322 :
1323 2025 : NULLIFY (cell)
1324 :
1325 : ! Subtract the vcom
1326 2025 : CALL compute_vcom(part, is_fixed, vcom)
1327 2025 : CALL subtract_vcom(part, is_fixed, vcom)
1328 : ! If requested and the system is not periodic, subtract the angular velocity
1329 2025 : CALL force_env_get(force_env, cell=cell)
1330 8100 : IF (SUM(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN
1331 4 : CALL compute_rcom(part, is_fixed, rcom)
1332 4 : CALL compute_vang(part, is_fixed, rcom, vang)
1333 4 : CALL subtract_vang(part, is_fixed, rcom, vang)
1334 : END IF
1335 : ! Rescale the velocities
1336 2025 : IF (simpar%do_thermal_region) THEN
1337 2 : CALL rescale_vel_region(part, md_env, simpar)
1338 : ELSE
1339 2023 : ekin = compute_ekin(part)
1340 2023 : CALL rescale_vel(part, simpar, ekin)
1341 : END IF
1342 2025 : END SUBROUTINE normalize_velocities
1343 :
1344 : ! **************************************************************************************************
1345 : !> \brief Computes Ekin, VCOM and Temp for particles
1346 : !> \param subsys ...
1347 : !> \param md_ener ...
1348 : !> \param vsubtract ...
1349 : !> \par History
1350 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1351 : ! **************************************************************************************************
1352 42 : SUBROUTINE reset_vcom(subsys, md_ener, vsubtract)
1353 : TYPE(cp_subsys_type), POINTER :: subsys
1354 : TYPE(md_ener_type), POINTER :: md_ener
1355 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vsubtract
1356 :
1357 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_vcom'
1358 :
1359 : INTEGER :: atom, handle, iatom, ikind, natom, &
1360 : shell_index
1361 42 : INTEGER, DIMENSION(:), POINTER :: atom_list
1362 : LOGICAL :: is_shell
1363 : REAL(KIND=dp) :: ekin_old, imass_c, imass_s, mass, v2
1364 : REAL(KIND=dp), DIMENSION(3) :: tmp, v
1365 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1366 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1367 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1368 : shell_particles
1369 : TYPE(shell_kind_type), POINTER :: shell
1370 :
1371 42 : NULLIFY (particles, atomic_kind, atomic_kinds, atom_list, shell)
1372 42 : CALL timeset(routineN, handle)
1373 :
1374 : CALL cp_subsys_get(subsys, &
1375 : atomic_kinds=atomic_kinds, &
1376 : particles=particles, &
1377 : shell_particles=shell_particles, &
1378 42 : core_particles=core_particles)
1379 :
1380 42 : ekin_old = md_ener%ekin
1381 : ! Possibly subtract a quantity from all velocities
1382 126 : DO ikind = 1, atomic_kinds%n_els
1383 84 : atomic_kind => atomic_kinds%els(ikind)
1384 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, &
1385 84 : natom=natom, mass=mass, shell_active=is_shell, shell=shell)
1386 126 : IF (is_shell) THEN
1387 336 : tmp = 0.5_dp*vsubtract*mass
1388 84 : imass_s = 1.0_dp/shell%mass_shell
1389 84 : imass_c = 1.0_dp/shell%mass_core
1390 3780 : DO iatom = 1, natom
1391 3696 : atom = atom_list(iatom)
1392 3696 : shell_index = particles%els(atom)%shell_index
1393 14784 : shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s
1394 14784 : core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c
1395 14868 : particles%els(atom)%v = particles%els(atom)%v - vsubtract
1396 : END DO
1397 : ELSE
1398 0 : DO iatom = 1, natom
1399 0 : atom = atom_list(iatom)
1400 0 : particles%els(atom)%v = particles%els(atom)%v - vsubtract
1401 : END DO
1402 : END IF
1403 : END DO
1404 : ! Compute Kinetic Energy and COM Velocity
1405 168 : md_ener%vcom = 0.0_dp
1406 42 : md_ener%total_mass = 0.0_dp
1407 42 : md_ener%ekin = 0.0_dp
1408 126 : DO ikind = 1, atomic_kinds%n_els
1409 84 : atomic_kind => atomic_kinds%els(ikind)
1410 84 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom)
1411 84 : v2 = 0.0_dp
1412 84 : v = 0.0_dp
1413 3780 : DO iatom = 1, natom
1414 3696 : atom = atom_list(iatom)
1415 14784 : v2 = v2 + SUM(particles%els(atom)%v**2)
1416 3696 : v(1) = v(1) + particles%els(atom)%v(1)
1417 3696 : v(2) = v(2) + particles%els(atom)%v(2)
1418 3780 : v(3) = v(3) + particles%els(atom)%v(3)
1419 : END DO
1420 84 : md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2
1421 84 : md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1)
1422 84 : md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2)
1423 84 : md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3)
1424 210 : md_ener%total_mass = md_ener%total_mass + REAL(natom, KIND=dp)*mass
1425 : END DO
1426 168 : md_ener%vcom = md_ener%vcom/md_ener%total_mass
1427 42 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1428 42 : IF (md_ener%nfree /= 0) THEN
1429 42 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1430 : END IF
1431 42 : CALL timestop(handle)
1432 :
1433 42 : END SUBROUTINE reset_vcom
1434 :
1435 : ! **************************************************************************************************
1436 : !> \brief Scale velocities to get the correct temperature
1437 : !> \param subsys ...
1438 : !> \param md_ener ...
1439 : !> \param temp_expected ...
1440 : !> \param temp_tol ...
1441 : !> \param iw ...
1442 : !> \par History
1443 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1444 : ! **************************************************************************************************
1445 12914 : SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
1446 : TYPE(cp_subsys_type), POINTER :: subsys
1447 : TYPE(md_ener_type), POINTER :: md_ener
1448 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1449 : INTEGER, INTENT(IN) :: iw
1450 :
1451 : REAL(KIND=dp) :: ekin_old, scale, temp_old
1452 :
1453 12914 : IF (ABS(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN
1454 2640 : scale = 0.0_dp
1455 2640 : IF (md_ener%temp_part > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_part)*kelvin)
1456 2640 : ekin_old = md_ener%ekin
1457 2640 : temp_old = md_ener%temp_part
1458 2640 : md_ener%ekin = 0.0_dp
1459 2640 : md_ener%temp_part = 0.0_dp
1460 10560 : md_ener%vcom = 0.0_dp
1461 2640 : md_ener%total_mass = 0.0_dp
1462 :
1463 2640 : CALL scale_velocity_low(subsys, scale, ireg=0, ekin=md_ener%ekin, vcom=md_ener%vcom)
1464 2640 : IF (md_ener%nfree /= 0) THEN
1465 2640 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1466 : END IF
1467 2640 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1468 2640 : IF (iw > 0) THEN
1469 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1470 1320 : 'MD_VEL| Temperature scaled to requested temperature'
1471 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1472 1320 : 'MD_VEL| Old temperature [K]', temp_old, &
1473 2640 : 'MD_VEL| New temperature [K]', md_ener%temp_part
1474 : END IF
1475 : END IF
1476 :
1477 12914 : END SUBROUTINE scale_velocity
1478 :
1479 : ! **************************************************************************************************
1480 : !> \brief Scale velocities of set of regions
1481 : !> \param md_env ...
1482 : !> \param subsys ...
1483 : !> \param md_ener ...
1484 : !> \param simpar ...
1485 : !> \param iw ...
1486 : !> \par author MI
1487 : ! **************************************************************************************************
1488 96 : SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1489 :
1490 : TYPE(md_environment_type), POINTER :: md_env
1491 : TYPE(cp_subsys_type), POINTER :: subsys
1492 : TYPE(md_ener_type), POINTER :: md_ener
1493 : TYPE(simpar_type), POINTER :: simpar
1494 : INTEGER, INTENT(IN) :: iw
1495 :
1496 : INTEGER :: ireg, nfree, nfree_done, nregions
1497 : REAL(KIND=dp) :: ekin, ekin_old, ekin_total_new, fscale, &
1498 : vcom(3), vcom_total(3)
1499 96 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp_new, temp_old
1500 : TYPE(particle_list_type), POINTER :: particles
1501 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1502 : TYPE(thermal_region_type), POINTER :: t_region
1503 : TYPE(thermal_regions_type), POINTER :: thermal_regions
1504 :
1505 96 : NULLIFY (particles, part, thermal_regions, t_region)
1506 96 : CALL cp_subsys_get(subsys, particles=particles)
1507 96 : part => particles%els
1508 96 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
1509 :
1510 96 : nregions = thermal_regions%nregions
1511 96 : nfree_done = 0
1512 96 : ekin_total_new = 0.0_dp
1513 96 : ekin_old = md_ener%ekin
1514 : vcom_total = 0.0_dp
1515 384 : ALLOCATE (temp_new(0:nregions), temp_old(0:nregions))
1516 96 : temp_new = 0.0_dp
1517 96 : temp_old = 0.0_dp
1518 : !loop regions
1519 256 : DO ireg = 1, nregions
1520 160 : NULLIFY (t_region)
1521 160 : t_region => thermal_regions%thermal_region(ireg)
1522 160 : nfree = 3*t_region%npart
1523 160 : ekin = compute_ekin(part, ireg)
1524 160 : IF (nfree > 0) t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1525 160 : temp_old(ireg) = t_region%temperature
1526 160 : IF (t_region%temp_tol > 0.0_dp .AND. &
1527 : ABS(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN
1528 2 : fscale = SQRT((t_region%temp_expected/t_region%temperature)*kelvin)
1529 2 : CALL scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1530 2 : t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1531 2 : temp_new(ireg) = t_region%temperature
1532 : END IF
1533 160 : nfree_done = nfree_done + nfree
1534 256 : ekin_total_new = ekin_total_new + ekin
1535 : END DO
1536 96 : nfree = simpar%nfree - nfree_done
1537 96 : ekin = compute_ekin(part, ireg=0)
1538 96 : IF (nfree > 0) thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1539 96 : temp_old(0) = thermal_regions%temp_reg0
1540 96 : IF (simpar%temp_tol > 0.0_dp .AND. nfree > 0) THEN
1541 0 : IF (ABS(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN
1542 0 : fscale = SQRT((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin)
1543 0 : CALL scale_velocity_low(subsys, fscale, 0, ekin, vcom)
1544 0 : thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1545 0 : temp_new(0) = thermal_regions%temp_reg0
1546 : END IF
1547 : END IF
1548 96 : ekin_total_new = ekin_total_new + ekin
1549 :
1550 96 : md_ener%ekin = ekin_total_new
1551 96 : IF (md_ener%nfree /= 0) THEN
1552 96 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1553 : END IF
1554 96 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1555 96 : IF (iw > 0) THEN
1556 176 : DO ireg = 0, nregions
1557 176 : IF (temp_new(ireg) > 0.0_dp) THEN
1558 : WRITE (UNIT=iw, FMT='(/,T2,A,I0,A)') &
1559 1 : 'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
1560 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1561 1 : 'MD_VEL| Old temperature [K]', temp_old(ireg), &
1562 2 : 'MD_VEL| New temperature [K]', temp_new(ireg)
1563 : END IF
1564 : END DO
1565 : END IF
1566 96 : DEALLOCATE (temp_new, temp_old)
1567 :
1568 96 : END SUBROUTINE scale_velocity_region
1569 :
1570 : ! **************************************************************************************************
1571 : !> \brief Scale velocities for a specific region
1572 : !> \param subsys ...
1573 : !> \param fscale ...
1574 : !> \param ireg ...
1575 : !> \param ekin ...
1576 : !> \param vcom ...
1577 : !> \par author MI
1578 : ! **************************************************************************************************
1579 2642 : SUBROUTINE scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1580 :
1581 : TYPE(cp_subsys_type), POINTER :: subsys
1582 : REAL(KIND=dp), INTENT(IN) :: fscale
1583 : INTEGER, INTENT(IN) :: ireg
1584 : REAL(KIND=dp), INTENT(OUT) :: ekin, vcom(3)
1585 :
1586 : INTEGER :: atom, iatom, ikind, my_ireg, natom, &
1587 : shell_index
1588 2642 : INTEGER, DIMENSION(:), POINTER :: atom_list
1589 : LOGICAL :: is_shell
1590 : REAL(KIND=dp) :: imass, mass, tmass, v2
1591 : REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1592 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1593 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1594 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1595 : shell_particles
1596 : TYPE(shell_kind_type), POINTER :: shell
1597 :
1598 2642 : NULLIFY (atomic_kinds, particles, shell_particles, core_particles, shell, atom_list)
1599 :
1600 2642 : my_ireg = ireg
1601 2642 : ekin = 0.0_dp
1602 2642 : tmass = 0.0_dp
1603 2642 : vcom = 0.0_dp
1604 :
1605 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
1606 2642 : shell_particles=shell_particles, core_particles=core_particles)
1607 :
1608 9082 : DO ikind = 1, atomic_kinds%n_els
1609 6440 : atomic_kind => atomic_kinds%els(ikind)
1610 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, &
1611 6440 : natom=natom, shell_active=is_shell, shell=shell)
1612 6440 : IF (is_shell) THEN
1613 124 : imass = 1.0_dp/mass
1614 124 : v2 = 0.0_dp
1615 124 : v = 0.0_dp
1616 5740 : DO iatom = 1, natom
1617 5616 : atom = atom_list(iatom)
1618 : !check region
1619 5616 : IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
1620 :
1621 21080 : particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1622 5270 : shell_index = particles%els(atom)%shell_index
1623 21080 : vs = shell_particles%els(shell_index)%v
1624 21080 : vc = core_particles%els(shell_index)%v
1625 5270 : tmp(1) = imass*(vs(1) - vc(1))
1626 5270 : tmp(2) = imass*(vs(2) - vc(2))
1627 5270 : tmp(3) = imass*(vs(3) - vc(3))
1628 :
1629 5270 : shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core
1630 5270 : shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core
1631 5270 : shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core
1632 :
1633 5270 : core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell
1634 5270 : core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell
1635 5270 : core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell
1636 :
1637 : ! kinetic energy and velocity of COM
1638 21080 : v2 = v2 + SUM(particles%els(atom)%v**2)
1639 5270 : v(1) = v(1) + particles%els(atom)%v(1)
1640 5270 : v(2) = v(2) + particles%els(atom)%v(2)
1641 5270 : v(3) = v(3) + particles%els(atom)%v(3)
1642 5740 : tmass = tmass + mass
1643 : END DO
1644 : ELSE
1645 6316 : v2 = 0.0_dp
1646 6316 : v = 0.0_dp
1647 22826 : DO iatom = 1, natom
1648 16510 : atom = atom_list(iatom)
1649 : !check region
1650 16510 : IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
1651 :
1652 66040 : particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1653 : ! kinetic energy and velocity of COM
1654 66040 : v2 = v2 + SUM(particles%els(atom)%v**2)
1655 16510 : v(1) = v(1) + particles%els(atom)%v(1)
1656 16510 : v(2) = v(2) + particles%els(atom)%v(2)
1657 16510 : v(3) = v(3) + particles%els(atom)%v(3)
1658 22826 : tmass = tmass + mass
1659 : END DO
1660 : END IF
1661 6440 : ekin = ekin + 0.5_dp*mass*v2
1662 6440 : vcom(1) = vcom(1) + mass*v(1)
1663 6440 : vcom(2) = vcom(2) + mass*v(2)
1664 15522 : vcom(3) = vcom(3) + mass*v(3)
1665 :
1666 : END DO
1667 10568 : vcom = vcom/tmass
1668 :
1669 2642 : END SUBROUTINE scale_velocity_low
1670 :
1671 : ! **************************************************************************************************
1672 : !> \brief Scale internal motion of CORE-SHELL model to the correct temperature
1673 : !> \param subsys ...
1674 : !> \param md_ener ...
1675 : !> \param temp_expected ...
1676 : !> \param temp_tol ...
1677 : !> \param iw ...
1678 : !> \par History
1679 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1680 : ! **************************************************************************************************
1681 1060 : SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
1682 : TYPE(cp_subsys_type), POINTER :: subsys
1683 : TYPE(md_ener_type), POINTER :: md_ener
1684 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1685 : INTEGER, INTENT(IN) :: iw
1686 :
1687 : INTEGER :: atom, iatom, ikind, natom, shell_index
1688 1060 : INTEGER, DIMENSION(:), POINTER :: atom_list
1689 : LOGICAL :: is_shell
1690 : REAL(KIND=dp) :: ekin_shell_old, fac_mass, mass, scale, &
1691 : temp_shell_old, v2
1692 : REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1693 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1694 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1695 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1696 : shell_particles
1697 : TYPE(shell_kind_type), POINTER :: shell
1698 :
1699 1060 : NULLIFY (atom_list, atomic_kinds, atomic_kind, core_particles, particles, shell_particles, shell)
1700 1060 : IF (ABS(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN
1701 80 : scale = 0.0_dp
1702 80 : IF (md_ener%temp_shell > EPSILON(0.0_dp)) scale = SQRT((temp_expected/md_ener%temp_shell)*kelvin)
1703 80 : ekin_shell_old = md_ener%ekin_shell
1704 80 : temp_shell_old = md_ener%temp_shell
1705 80 : md_ener%ekin_shell = 0.0_dp
1706 80 : md_ener%temp_shell = 0.0_dp
1707 :
1708 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles, &
1709 80 : core_particles=core_particles)
1710 :
1711 240 : DO ikind = 1, atomic_kinds%n_els
1712 160 : atomic_kind => atomic_kinds%els(ikind)
1713 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom, &
1714 160 : shell_active=is_shell, shell=shell)
1715 240 : IF (is_shell) THEN
1716 160 : fac_mass = 1.0_dp/mass
1717 160 : v2 = 0.0_dp
1718 776 : DO iatom = 1, natom
1719 616 : atom = atom_list(iatom)
1720 616 : shell_index = particles%els(atom)%shell_index
1721 2464 : vs = shell_particles%els(shell_index)%v
1722 2464 : vc = core_particles%els(shell_index)%v
1723 2464 : v = particles%els(atom)%v
1724 616 : tmp(1) = fac_mass*(vc(1) - vs(1))
1725 616 : tmp(2) = fac_mass*(vc(2) - vs(2))
1726 616 : tmp(3) = fac_mass*(vc(3) - vs(3))
1727 :
1728 616 : shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1)
1729 616 : shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2)
1730 616 : shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3)
1731 :
1732 616 : core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1)
1733 616 : core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2)
1734 616 : core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3)
1735 :
1736 2464 : vs = shell_particles%els(shell_index)%v
1737 2464 : vc = core_particles%els(shell_index)%v
1738 616 : tmp(1) = vc(1) - vs(1)
1739 616 : tmp(2) = vc(2) - vs(2)
1740 616 : tmp(3) = vc(3) - vs(3)
1741 2624 : v2 = v2 + SUM(tmp**2)
1742 : END DO
1743 160 : md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2
1744 : END IF
1745 : END DO
1746 80 : IF (md_ener%nfree_shell > 0) THEN
1747 80 : md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)*kelvin
1748 : END IF
1749 80 : md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell
1750 80 : IF (iw > 0) THEN
1751 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1752 40 : 'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
1753 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1754 40 : 'MD_VEL| Old temperature [K]', temp_shell_old, &
1755 80 : 'MD_VEL| New temperature [K]', md_ener%temp_shell
1756 : END IF
1757 : END IF
1758 :
1759 1060 : END SUBROUTINE scale_velocity_internal
1760 :
1761 : ! **************************************************************************************************
1762 : !> \brief Scale barostat velocities to get the desired temperature
1763 : !> \param md_env ...
1764 : !> \param md_ener ...
1765 : !> \param temp_expected ...
1766 : !> \param temp_tol ...
1767 : !> \param iw ...
1768 : !> \par History
1769 : !> MI 02.2008
1770 : ! **************************************************************************************************
1771 40 : SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
1772 : TYPE(md_environment_type), POINTER :: md_env
1773 : TYPE(md_ener_type), POINTER :: md_ener
1774 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1775 : INTEGER, INTENT(IN) :: iw
1776 :
1777 : INTEGER :: i, j, nfree
1778 : REAL(KIND=dp) :: ekin_old, scale, temp_old
1779 40 : TYPE(npt_info_type), POINTER :: npt(:, :)
1780 : TYPE(simpar_type), POINTER :: simpar
1781 :
1782 40 : NULLIFY (npt, simpar)
1783 40 : CALL get_md_env(md_env, simpar=simpar, npt=npt)
1784 40 : IF (ABS(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN
1785 2 : scale = 0.0_dp
1786 2 : IF (md_ener%temp_baro > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_baro)*kelvin)
1787 2 : ekin_old = md_ener%baro_kin
1788 2 : temp_old = md_ener%temp_baro
1789 2 : md_ener%baro_kin = 0.0_dp
1790 2 : md_ener%temp_baro = 0.0_dp
1791 : IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
1792 2 : .OR. simpar%ensemble == npt_ia_ensemble) THEN
1793 0 : npt(1, 1)%v = npt(1, 1)%v*scale
1794 0 : md_ener%baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
1795 : ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
1796 2 : md_ener%baro_kin = 0.0_dp
1797 8 : DO i = 1, 3
1798 26 : DO j = 1, 3
1799 18 : npt(i, j)%v = npt(i, j)%v*scale
1800 24 : md_ener%baro_kin = md_ener%baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
1801 : END DO
1802 : END DO
1803 : END IF
1804 2 : nfree = SIZE(npt, 1)*SIZE(npt, 2)
1805 2 : md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/REAL(nfree, dp)*kelvin
1806 2 : IF (iw > 0) THEN
1807 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1808 1 : 'MD_VEL| Temperature of barostat motion scaled to requested temperature'
1809 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1810 1 : 'MD_VEL| Old temperature [K]', temp_old, &
1811 2 : 'MD_VEL| New temperature [K]', md_ener%temp_baro
1812 : END IF
1813 : END IF
1814 :
1815 40 : END SUBROUTINE scale_velocity_baro
1816 :
1817 : ! **************************************************************************************************
1818 : !> \brief Perform all temperature manipulations during a QS MD run.
1819 : !> \param simpar ...
1820 : !> \param md_env ...
1821 : !> \param md_ener ...
1822 : !> \param force_env ...
1823 : !> \param logger ...
1824 : !> \par History
1825 : !> Creation (15.09.2003,MK)
1826 : !> adapted to force_env (05.10.2003,fawzi)
1827 : !> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1828 : ! **************************************************************************************************
1829 40248 : SUBROUTINE temperature_control(simpar, md_env, md_ener, force_env, logger)
1830 :
1831 : TYPE(simpar_type), POINTER :: simpar
1832 : TYPE(md_environment_type), POINTER :: md_env
1833 : TYPE(md_ener_type), POINTER :: md_ener
1834 : TYPE(force_env_type), POINTER :: force_env
1835 : TYPE(cp_logger_type), POINTER :: logger
1836 :
1837 : CHARACTER(LEN=*), PARAMETER :: routineN = 'temperature_control'
1838 :
1839 : INTEGER :: handle, iw
1840 : TYPE(cp_subsys_type), POINTER :: subsys
1841 : TYPE(mp_para_env_type), POINTER :: para_env
1842 :
1843 40248 : CALL timeset(routineN, handle)
1844 40248 : NULLIFY (subsys, para_env)
1845 40248 : CPASSERT(ASSOCIATED(simpar))
1846 40248 : CPASSERT(ASSOCIATED(md_ener))
1847 40248 : CPASSERT(ASSOCIATED(force_env))
1848 40248 : CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
1849 : iw = cp_print_key_unit_nr(logger, force_env%root_section, "MOTION%MD%PRINT%PROGRAM_RUN_INFO", &
1850 40248 : extension=".mdLog")
1851 :
1852 : ! Control the particle motion
1853 40248 : IF (simpar%do_thermal_region) THEN
1854 96 : CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1855 : ELSE
1856 40152 : IF (simpar%temp_tol > 0.0_dp) THEN
1857 12870 : CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw)
1858 : END IF
1859 : END IF
1860 : ! Control the internal core-shell motion
1861 40248 : IF (simpar%temp_sh_tol > 0.0_dp) THEN
1862 1060 : CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw)
1863 : END IF
1864 : ! Control cell motion
1865 42768 : SELECT CASE (simpar%ensemble)
1866 : CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
1867 : npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_ia_ensemble)
1868 40248 : IF (simpar%temp_baro_tol > 0.0_dp) THEN
1869 40 : CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw)
1870 : END IF
1871 : END SELECT
1872 :
1873 : CALL cp_print_key_finished_output(iw, logger, force_env%root_section, &
1874 40248 : "MOTION%MD%PRINT%PROGRAM_RUN_INFO")
1875 40248 : CALL timestop(handle)
1876 40248 : END SUBROUTINE temperature_control
1877 :
1878 : ! **************************************************************************************************
1879 : !> \brief Set to 0 the velocity of the COM along MD runs, if required.
1880 : !> \param md_ener ...
1881 : !> \param force_env ...
1882 : !> \param md_section ...
1883 : !> \param logger ...
1884 : !> \par History
1885 : !> Creation (29.04.2007,MI)
1886 : !> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1887 : ! **************************************************************************************************
1888 80496 : SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
1889 :
1890 : TYPE(md_ener_type), POINTER :: md_ener
1891 : TYPE(force_env_type), POINTER :: force_env
1892 : TYPE(section_vals_type), POINTER :: md_section
1893 : TYPE(cp_logger_type), POINTER :: logger
1894 :
1895 : CHARACTER(LEN=*), PARAMETER :: routineN = 'comvel_control'
1896 :
1897 : INTEGER :: handle, iw
1898 : LOGICAL :: explicit
1899 : REAL(KIND=dp) :: comvel_tol, temp_old, vel_com
1900 : REAL(KIND=dp), DIMENSION(3) :: vcom_old
1901 : TYPE(cp_subsys_type), POINTER :: subsys
1902 :
1903 40248 : CALL timeset(routineN, handle)
1904 40248 : NULLIFY (subsys)
1905 40248 : CPASSERT(ASSOCIATED(force_env))
1906 40248 : CALL force_env_get(force_env, subsys=subsys)
1907 :
1908 : ! Print COMVEL and COM Position
1909 40248 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
1910 40248 : IF (iw > 0) THEN
1911 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1912 7785 : "MD_VEL| Centre of mass motion (COM)"
1913 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1914 31140 : "MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
1915 : END IF
1916 40248 : CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")
1917 :
1918 : ! If requested rescale COMVEL
1919 40248 : CALL section_vals_val_get(md_section, "COMVEL_TOL", explicit=explicit)
1920 40248 : IF (explicit) THEN
1921 826 : CALL section_vals_val_get(md_section, "COMVEL_TOL", r_val=comvel_tol)
1922 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1923 826 : extension=".mdLog")
1924 826 : vel_com = SQRT(md_ener%vcom(1)**2 + md_ener%vcom(2)**2 + md_ener%vcom(3)**2)
1925 :
1926 : ! Subtract the velocity of the COM, if requested
1927 826 : IF (vel_com > comvel_tol) THEN
1928 42 : temp_old = md_ener%temp_part/kelvin
1929 168 : vcom_old = md_ener%vcom
1930 42 : CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
1931 42 : CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
1932 42 : IF (iw > 0) THEN
1933 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1934 21 : "MD_VEL| Old VCOM [a.u.]", vcom_old(1:3)
1935 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1936 84 : "MD_VEL| New VCOM [a.u.]", md_ener%vcom(1:3)
1937 : END IF
1938 : END IF
1939 : CALL cp_print_key_finished_output(iw, logger, md_section, &
1940 826 : "PRINT%PROGRAM_RUN_INFO")
1941 : END IF
1942 :
1943 40248 : CALL timestop(handle)
1944 40248 : END SUBROUTINE comvel_control
1945 :
1946 : ! **************************************************************************************************
1947 : !> \brief Set to 0 the angular velocity along MD runs, if required.
1948 : !> \param md_ener ...
1949 : !> \param force_env ...
1950 : !> \param md_section ...
1951 : !> \param logger ...
1952 : !> \par History
1953 : !> Creation (10.2009) Teodoro Laino [tlaino]
1954 : ! **************************************************************************************************
1955 40248 : SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
1956 :
1957 : TYPE(md_ener_type), POINTER :: md_ener
1958 : TYPE(force_env_type), POINTER :: force_env
1959 : TYPE(section_vals_type), POINTER :: md_section
1960 : TYPE(cp_logger_type), POINTER :: logger
1961 :
1962 : CHARACTER(LEN=*), PARAMETER :: routineN = 'angvel_control'
1963 :
1964 : INTEGER :: handle, ifixd, imolecule_kind, iw, natoms
1965 40248 : INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
1966 : LOGICAL :: explicit
1967 : REAL(KIND=dp) :: angvel_tol, rcom(3), temp_old, vang(3), &
1968 : vang_new(3)
1969 : TYPE(cell_type), POINTER :: cell
1970 : TYPE(cp_subsys_type), POINTER :: subsys
1971 40248 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
1972 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1973 40248 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1974 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1975 : TYPE(particle_list_type), POINTER :: particles
1976 :
1977 40248 : CALL timeset(routineN, handle)
1978 : ! If requested rescale ANGVEL
1979 40248 : CALL section_vals_val_get(md_section, "ANGVEL_TOL", explicit=explicit)
1980 40248 : IF (explicit) THEN
1981 40 : NULLIFY (subsys, cell)
1982 40 : CPASSERT(ASSOCIATED(force_env))
1983 40 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
1984 :
1985 160 : IF (SUM(cell%perd(1:3)) == 0) THEN
1986 40 : CALL section_vals_val_get(md_section, "ANGVEL_TOL", r_val=angvel_tol)
1987 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1988 40 : extension=".mdLog")
1989 :
1990 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, &
1991 40 : particles=particles)
1992 :
1993 40 : natoms = SIZE(particles%els)
1994 : ! Build a list of all fixed atoms (if any)
1995 120 : ALLOCATE (is_fixed(natoms))
1996 :
1997 40 : is_fixed = use_perd_none
1998 40 : molecule_kind_set => molecule_kinds%els
1999 600 : DO imolecule_kind = 1, molecule_kinds%n_els
2000 560 : molecule_kind => molecule_kind_set(imolecule_kind)
2001 560 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2002 600 : IF (ASSOCIATED(fixd_list)) THEN
2003 0 : DO ifixd = 1, SIZE(fixd_list)
2004 0 : IF (.NOT. fixd_list(ifixd)%restraint%active) &
2005 0 : is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2006 : END DO
2007 : END IF
2008 : END DO
2009 :
2010 : ! If requested and the system is not periodic, subtract the angular velocity
2011 40 : CALL compute_rcom(particles%els, is_fixed, rcom)
2012 40 : CALL compute_vang(particles%els, is_fixed, rcom, vang)
2013 : ! SQRT(DOT_PRODUCT(vang,vang))>angvel_tol
2014 160 : IF (DOT_PRODUCT(vang, vang) > (angvel_tol*angvel_tol)) THEN
2015 2 : CALL subtract_vang(particles%els, is_fixed, rcom, vang)
2016 :
2017 : ! Rescale velocities after removal
2018 2 : temp_old = md_ener%temp_part/kelvin
2019 2 : CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
2020 2 : CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
2021 2 : IF (iw > 0) THEN
2022 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
2023 1 : 'MD_VEL| Old VANG [a.u.]', vang(1:3)
2024 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
2025 1 : 'MD_VEL| New VANG [a.u.]', vang_new(1:3)
2026 : END IF
2027 : END IF
2028 :
2029 40 : DEALLOCATE (is_fixed)
2030 :
2031 : CALL cp_print_key_finished_output(iw, logger, md_section, &
2032 80 : "PRINT%PROGRAM_RUN_INFO")
2033 : END IF
2034 : END IF
2035 :
2036 40248 : CALL timestop(handle)
2037 80496 : END SUBROUTINE angvel_control
2038 :
2039 : ! **************************************************************************************************
2040 : !> \brief Initialize Velocities for MD runs
2041 : !> \param force_env ...
2042 : !> \param simpar ...
2043 : !> \param globenv ...
2044 : !> \param md_env ...
2045 : !> \param md_section ...
2046 : !> \param constraint_section ...
2047 : !> \param write_binary_restart_file ...
2048 : !> \par History
2049 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
2050 : ! **************************************************************************************************
2051 3534 : SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, &
2052 : constraint_section, write_binary_restart_file)
2053 :
2054 : TYPE(force_env_type), POINTER :: force_env
2055 : TYPE(simpar_type), POINTER :: simpar
2056 : TYPE(global_environment_type), POINTER :: globenv
2057 : TYPE(md_environment_type), POINTER :: md_env
2058 : TYPE(section_vals_type), POINTER :: md_section, constraint_section
2059 : LOGICAL, INTENT(IN) :: write_binary_restart_file
2060 :
2061 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_velocities'
2062 :
2063 : INTEGER :: handle, nconstraint, nconstraint_fixd
2064 : LOGICAL :: apply_cns0, shell_adiabatic, &
2065 : shell_present
2066 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2067 : TYPE(cell_type), POINTER :: cell
2068 : TYPE(cp_subsys_type), POINTER :: subsys
2069 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2070 : TYPE(mp_para_env_type), POINTER :: para_env
2071 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
2072 : shell_particles
2073 1767 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
2074 1767 : shell_particle_set
2075 : TYPE(section_vals_type), POINTER :: force_env_section, print_section, &
2076 : subsys_section
2077 :
2078 1767 : CALL timeset(routineN, handle)
2079 :
2080 1767 : NULLIFY (atomic_kinds, cell, para_env, subsys, molecule_kinds, core_particles, particles)
2081 1767 : NULLIFY (shell_particles, core_particle_set, particle_set, shell_particle_set)
2082 1767 : NULLIFY (force_env_section, print_section, subsys_section)
2083 :
2084 1767 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
2085 1767 : apply_cns0 = .FALSE.
2086 1767 : IF (simpar%constraint) THEN
2087 310 : CALL section_vals_val_get(constraint_section, "CONSTRAINT_INIT", l_val=apply_cns0)
2088 : END IF
2089 : ! Always initialize velocities and possibly restart them
2090 : CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env, &
2091 1767 : force_env_section=force_env_section)
2092 1767 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2093 :
2094 : CALL cp_subsys_get(subsys, &
2095 : atomic_kinds=atomic_kinds, &
2096 : core_particles=core_particles, &
2097 : molecule_kinds=molecule_kinds, &
2098 : particles=particles, &
2099 1767 : shell_particles=shell_particles)
2100 :
2101 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
2102 : shell_present=shell_present, &
2103 1767 : shell_adiabatic=shell_adiabatic)
2104 :
2105 1767 : NULLIFY (core_particle_set)
2106 : NULLIFY (particle_set)
2107 1767 : NULLIFY (shell_particle_set)
2108 1767 : particle_set => particles%els
2109 :
2110 1767 : IF (shell_present .AND. shell_adiabatic) THEN
2111 : ! Constraints are not yet implemented for core-shell models generally
2112 : CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els, &
2113 : nconstraint=nconstraint, &
2114 132 : nconstraint_fixd=nconstraint_fixd)
2115 132 : IF (nconstraint - nconstraint_fixd /= 0) &
2116 0 : CPABORT("Only the fixed atom constraint is implemented for core-shell models")
2117 : !MK CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,failure)
2118 132 : CPASSERT(ASSOCIATED(shell_particles))
2119 132 : CPASSERT(ASSOCIATED(core_particles))
2120 132 : shell_particle_set => shell_particles%els
2121 132 : core_particle_set => core_particles%els
2122 : END IF
2123 :
2124 : CALL initialize_velocities(simpar, &
2125 : particle_set, &
2126 : molecule_kinds=molecule_kinds, &
2127 : force_env=force_env, &
2128 : globenv=globenv, &
2129 : md_env=md_env, &
2130 : label="Velocities initialization", &
2131 : print_section=print_section, &
2132 : subsys_section=subsys_section, &
2133 : shell_present=(shell_present .AND. shell_adiabatic), &
2134 : shell_part=shell_particle_set, &
2135 : core_part=core_particle_set, &
2136 : force_rescaling=.FALSE., &
2137 : para_env=para_env, &
2138 3402 : write_binary_restart_file=write_binary_restart_file)
2139 :
2140 : ! Apply constraints if required and rescale velocities..
2141 1767 : IF (simpar%ensemble /= reftraj_ensemble) THEN
2142 1729 : IF (apply_cns0) THEN
2143 24 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
2144 : CALL force_env_shake(force_env, &
2145 : shake_tol=simpar%shake_tol, &
2146 : log_unit=simpar%info_constraint, &
2147 : lagrange_mult=simpar%lagrange_multipliers, &
2148 : dump_lm=simpar%dump_lm, &
2149 24 : compold=.TRUE.)
2150 : CALL force_env_rattle(force_env, shake_tol=simpar%shake_tol, &
2151 : log_unit=simpar%info_constraint, lagrange_mult=simpar%lagrange_multipliers, &
2152 24 : dump_lm=simpar%dump_lm, reset=.TRUE.)
2153 24 : IF (simpar%do_respa) THEN
2154 : CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
2155 0 : calc_force=.TRUE.)
2156 : CALL force_env_shake(force_env%sub_force_env(1)%force_env, &
2157 : shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2158 0 : lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, compold=.TRUE.)
2159 : CALL force_env_rattle(force_env%sub_force_env(1)%force_env, &
2160 : shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2161 0 : lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, reset=.TRUE.)
2162 : END IF
2163 : ! Reinitialize velocities rescaling properly after rattle
2164 24 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2165 24 : CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
2166 : CALL initialize_velocities(simpar, &
2167 : particle_set, &
2168 : molecule_kinds=molecule_kinds, &
2169 : force_env=force_env, &
2170 : globenv=globenv, &
2171 : md_env=md_env, &
2172 : label="Re-Initializing velocities after applying constraints", &
2173 : print_section=print_section, &
2174 : subsys_section=subsys_section, &
2175 : shell_present=(shell_present .AND. shell_adiabatic), &
2176 : shell_part=shell_particle_set, &
2177 : core_part=core_particle_set, &
2178 : force_rescaling=.TRUE., &
2179 : para_env=para_env, &
2180 48 : write_binary_restart_file=write_binary_restart_file)
2181 : END IF
2182 : END IF
2183 :
2184 : ! Perform setup for a cascade run
2185 1767 : CALL initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2186 :
2187 1767 : CALL timestop(handle)
2188 :
2189 1767 : END SUBROUTINE setup_velocities
2190 :
2191 : ! **************************************************************************************************
2192 : !> \brief Perform the initialization for a cascade run
2193 : !> \param simpar ...
2194 : !> \param particle_set ...
2195 : !> \param molecule_kinds ...
2196 : !> \param md_section ...
2197 : !> \date 05.02.2012
2198 : !> \author Matthias Krack (MK)
2199 : !> \version 1.0
2200 : ! **************************************************************************************************
2201 1767 : SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2202 :
2203 : TYPE(simpar_type), POINTER :: simpar
2204 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2205 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2206 : TYPE(section_vals_type), POINTER :: md_section
2207 :
2208 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_cascade'
2209 :
2210 : CHARACTER(len=2*default_string_length) :: line
2211 : INTEGER :: handle, iatom, ifixd, imolecule_kind, &
2212 : iparticle, iw, natom, nparticle
2213 1767 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_index, is_fixed
2214 : LOGICAL :: init_cascade, is_ok, no_read_error
2215 : REAL(KIND=dp) :: ecom, ekin, energy, norm, temp, &
2216 : temperature
2217 1767 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: matom, weight
2218 1767 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vatom
2219 : REAL(KIND=dp), DIMENSION(3) :: vcom
2220 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2221 : TYPE(cp_logger_type), POINTER :: logger
2222 : TYPE(cp_sll_val_type), POINTER :: atom_list
2223 1767 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
2224 1767 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2225 : TYPE(molecule_kind_type), POINTER :: molecule_kind
2226 : TYPE(section_vals_type), POINTER :: atom_list_section, cascade_section, &
2227 : print_section
2228 : TYPE(val_type), POINTER :: val
2229 :
2230 1767 : CALL timeset(routineN, handle)
2231 :
2232 1767 : NULLIFY (atom_list)
2233 1767 : NULLIFY (atom_list_section)
2234 1767 : NULLIFY (atomic_kind)
2235 1767 : NULLIFY (cascade_section)
2236 1767 : NULLIFY (fixd_list)
2237 1767 : NULLIFY (molecule_kind)
2238 1767 : NULLIFY (molecule_kind_set)
2239 1767 : NULLIFY (logger)
2240 1767 : NULLIFY (val)
2241 :
2242 1767 : logger => cp_get_default_logger()
2243 1767 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
2244 1767 : iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
2245 :
2246 1767 : cascade_section => section_vals_get_subs_vals(md_section, "CASCADE")
2247 1767 : CALL section_vals_val_get(cascade_section, "_SECTION_PARAMETERS_", l_val=init_cascade)
2248 :
2249 1767 : nparticle = SIZE(particle_set)
2250 :
2251 1767 : IF (init_cascade) THEN
2252 :
2253 2 : CALL section_vals_val_get(cascade_section, "ENERGY", r_val=energy)
2254 2 : IF (energy < 0.0_dp) &
2255 0 : CPABORT("Error occurred reading &CASCADE section: Negative energy found")
2256 :
2257 2 : IF (iw > 0) THEN
2258 1 : ekin = cp_unit_from_cp2k(energy, "keV")
2259 : WRITE (UNIT=iw, FMT="(/,T2,A,T61,F20.6)") &
2260 1 : "CASCADE| Energy [keV]", ekin
2261 : WRITE (UNIT=iw, FMT="(T2,A)") &
2262 1 : "CASCADE|"
2263 : END IF
2264 :
2265 : ! Read the atomic velocities given in the input file
2266 2 : atom_list_section => section_vals_get_subs_vals(cascade_section, "ATOM_LIST")
2267 2 : CALL section_vals_val_get(atom_list_section, "_DEFAULT_KEYWORD_", n_rep_val=natom)
2268 2 : CALL section_vals_list_get(atom_list_section, "_DEFAULT_KEYWORD_", list=atom_list)
2269 2 : IF (natom <= 0) &
2270 0 : CPABORT("Error occurred reading &CASCADE section: No atom list found")
2271 :
2272 2 : IF (iw > 0) THEN
2273 : WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
2274 1 : "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "weight"
2275 : END IF
2276 :
2277 6 : ALLOCATE (atom_index(natom))
2278 6 : ALLOCATE (matom(natom))
2279 6 : ALLOCATE (vatom(3, natom))
2280 4 : ALLOCATE (weight(natom))
2281 :
2282 8 : DO iatom = 1, natom
2283 6 : is_ok = cp_sll_val_next(atom_list, val)
2284 6 : CALL val_get(val, c_val=line)
2285 : ! Read atomic index, velocity vector, and weight
2286 6 : no_read_error = .FALSE.
2287 6 : READ (UNIT=line, FMT=*, ERR=999) atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2288 : no_read_error = .TRUE.
2289 : 999 IF (.NOT. no_read_error) &
2290 0 : CPABORT("Error occurred reading &CASCADE section. Last line read <"//TRIM(line)//">")
2291 6 : IF ((atom_index(iatom) <= 0) .OR. ((atom_index(iatom) > nparticle))) &
2292 0 : CPABORT("Error occurred reading &CASCADE section: Invalid atom index found")
2293 6 : IF (weight(iatom) < 0.0_dp) &
2294 0 : CPABORT("Error occurred reading &CASCADE section: Negative weight found")
2295 8 : IF (iw > 0) THEN
2296 : WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
2297 3 : "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2298 : END IF
2299 : END DO
2300 :
2301 : ! Normalise velocities and weights
2302 : norm = 0.0_dp
2303 8 : DO iatom = 1, natom
2304 6 : iparticle = atom_index(iatom)
2305 6 : IF (particle_set(iparticle)%shell_index /= 0) THEN
2306 0 : CPWARN("Warning: The primary knock-on atom is a core-shell atom")
2307 : END IF
2308 6 : atomic_kind => particle_set(iparticle)%atomic_kind
2309 6 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
2310 8 : norm = norm + matom(iatom)*weight(iatom)
2311 : END DO
2312 8 : weight(:) = matom(:)*weight(:)*energy/norm
2313 8 : DO iatom = 1, natom
2314 24 : norm = SQRT(DOT_PRODUCT(vatom(1:3, iatom), vatom(1:3, iatom)))
2315 26 : vatom(1:3, iatom) = vatom(1:3, iatom)/norm
2316 : END DO
2317 :
2318 2 : IF (iw > 0) THEN
2319 : WRITE (UNIT=iw, FMT="(T2,A)") &
2320 1 : "CASCADE|", &
2321 1 : "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
2322 2 : "CASCADE|"
2323 : WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
2324 1 : "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
2325 4 : DO iatom = 1, natom
2326 3 : ekin = cp_unit_from_cp2k(weight(iatom), "keV")
2327 : WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
2328 4 : "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
2329 : END DO
2330 : END IF
2331 :
2332 : ! Apply velocity modifications
2333 8 : DO iatom = 1, natom
2334 6 : iparticle = atom_index(iatom)
2335 : particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
2336 26 : SQRT(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
2337 : END DO
2338 :
2339 2 : DEALLOCATE (atom_index)
2340 2 : DEALLOCATE (matom)
2341 2 : DEALLOCATE (vatom)
2342 2 : DEALLOCATE (weight)
2343 :
2344 6 : IF (iw > 0) THEN
2345 : ! Build a list of all fixed atoms (if any)
2346 3 : ALLOCATE (is_fixed(nparticle))
2347 1 : is_fixed = use_perd_none
2348 1 : molecule_kind_set => molecule_kinds%els
2349 2 : DO imolecule_kind = 1, molecule_kinds%n_els
2350 1 : molecule_kind => molecule_kind_set(imolecule_kind)
2351 1 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2352 2 : IF (ASSOCIATED(fixd_list)) THEN
2353 0 : DO ifixd = 1, SIZE(fixd_list)
2354 0 : IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2355 : END DO
2356 : END IF
2357 : END DO
2358 : ! Compute vcom, ecom and ekin for printout
2359 1 : CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
2360 1 : ekin = compute_ekin(particle_set) - ecom
2361 1 : IF (simpar%nfree == 0) THEN
2362 0 : CPASSERT(ekin == 0.0_dp)
2363 0 : temp = 0.0_dp
2364 : ELSE
2365 1 : temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
2366 : END IF
2367 1 : temperature = cp_unit_from_cp2k(temp, "K")
2368 : WRITE (UNIT=iw, FMT="(T2,A)") &
2369 1 : "CASCADE|"
2370 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.6)") &
2371 1 : "CASCADE| Temperature after cascade initialization [K]", temperature
2372 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,ES16.8))") &
2373 1 : "CASCADE| COM velocity", vcom(1:3)
2374 : !MK ! compute and log rcom and vang if not periodic
2375 : !MK CALL force_env_get(force_env,cell=cell)
2376 : !MK IF (SUM(cell%perd(1:3)) == 0) THEN
2377 : !MK CALL compute_rcom(particle_set,is_fixed,rcom)
2378 : !MK CALL compute_vang(particle_set,is_fixed,rcom,vang)
2379 : !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
2380 : !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
2381 : !MK END IF
2382 2 : DEALLOCATE (is_fixed)
2383 : END IF
2384 :
2385 : END IF
2386 :
2387 1767 : CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
2388 :
2389 1767 : CALL timestop(handle)
2390 :
2391 3534 : END SUBROUTINE initialize_cascade
2392 :
2393 : END MODULE md_vel_utils
|
