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 3179 : 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 3179 : vcom = 0.0_dp
152 3179 : denom = 0.0_dp
153 747154 : DO i = 1, SIZE(part)
154 743975 : atomic_kind => part(i)%atomic_kind
155 743975 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
156 747154 : IF (mass /= 0.0) THEN
157 1480582 : 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 737933 : vcom(1) = vcom(1) + part(i)%v(1)*mass
160 737933 : vcom(2) = vcom(2) + part(i)%v(2)*mass
161 737933 : vcom(3) = vcom(3) + part(i)%v(3)*mass
162 742649 : denom = denom + mass
163 : END SELECT
164 : END IF
165 : END DO
166 12716 : vcom = vcom/denom
167 3179 : IF (PRESENT(ecom)) THEN
168 4364 : ecom = 0.5_dp*denom*SUM(vcom*vcom)
169 : END IF
170 :
171 3179 : 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 3445 : 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 3445 : NULLIFY (atomic_kind)
222 3445 : ekin = 0.0_dp
223 3445 : 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 746960 : DO i = 1, SIZE(part)
233 743783 : atomic_kind => part(i)%atomic_kind
234 743783 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
235 2978309 : ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
236 : END DO
237 : END IF
238 :
239 3445 : 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 2116 : 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 2116 : 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 2110 : ELSE IF (PRESENT(nfree)) THEN
273 0 : my_ireg = 0
274 0 : my_nfree = nfree
275 0 : my_temp = simpar%temp_ext
276 : ELSE
277 2110 : my_ireg = 0
278 2110 : my_nfree = simpar%nfree
279 2110 : my_temp = simpar%temp_ext
280 : END IF
281 2116 : IF (my_nfree /= 0) THEN
282 2104 : 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 2116 : ekin = ekin*factor
290 2116 : factor = SQRT(factor)
291 2116 : 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 445292 : DO i = 1, SIZE(part)
297 1774838 : part(i)%v(:) = factor*part(i)%v(:)
298 : END DO
299 2110 : IF (PRESENT(vcom)) THEN
300 96 : vcom = factor*vcom
301 : END IF
302 : END IF
303 :
304 2116 : 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 2088 : 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 443592 : DO i = 1, SIZE(part)
370 2088 : 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 1750518 : part(i)%v(:) = part(i)%v(:) - vcom(:)
388 : END SELECT
389 : END DO
390 2088 : 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 1794 : 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 1794 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
570 1794 : 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 1794 : CALL timeset(routineN, handle)
575 :
576 : ! Initializing parameters
577 1794 : natoms = SIZE(part)
578 1794 : NULLIFY (atomic_kind, fixd_list, logger, molecule_kind)
579 1794 : NULLIFY (molecule_kind_set)
580 :
581 : ! Logging
582 1794 : logger => cp_get_default_logger()
583 1794 : 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 5382 : ALLOCATE (is_fixed(natoms))
587 :
588 1794 : is_fixed = use_perd_none
589 1794 : molecule_kind_set => molecule_kinds%els
590 66614 : DO imolecule_kind = 1, molecule_kinds%n_els
591 64820 : molecule_kind => molecule_kind_set(imolecule_kind)
592 64820 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
593 66614 : 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 1794 : 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 1794 : shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
614 1794 : IF (.NOT. success) THEN
615 3058 : SELECT CASE (simpar%initialization_method)
616 : CASE (md_init_default)
617 : CALL generate_velocities(simpar, part, force_env, globenv, md_env, shell_present, &
618 1528 : 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 1532 : CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
635 : END SELECT
636 : END IF
637 :
638 1794 : 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 1794 : DEALLOCATE (is_fixed)
669 1794 : CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
670 1794 : CALL timestop(handle)
671 :
672 3588 : 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 12558 : SUBROUTINE read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
691 1794 : 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 1794 : 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 1794 : success = .FALSE.
722 1794 : natoms = SIZE(part)
723 : atomvel_read = .FALSE.
724 : corevel_read = .FALSE.
725 : shellvel_read = .FALSE.
726 1794 : NULLIFY (vel, atomic_kind, atom_list, core_list, shell_list)
727 1794 : NULLIFY (atomvel_section, shellvel_section, corevel_section)
728 1794 : NULLIFY (cell, shell, thermal_regions, val)
729 1794 : CALL force_env_get(force_env, cell=cell)
730 :
731 : ! Core-Shell Model
732 1794 : nshell = 0
733 1794 : 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 1794 : atomvel_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
740 1794 : shellvel_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
741 1794 : corevel_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
742 :
743 : ! Read or initialize the particle velocities
744 1794 : CALL section_vals_get(atomvel_section, explicit=atomvel_explicit)
745 1794 : CALL section_vals_get(shellvel_section, explicit=shellvel_explicit)
746 1794 : CALL section_vals_get(corevel_section, explicit=corevel_explicit)
747 1794 : CPASSERT(shellvel_explicit .EQV. corevel_explicit)
748 :
749 : CALL read_binary_velocities("", part, force_env%root_section, para_env, &
750 1794 : subsys_section, atomvel_read, cell)
751 : CALL read_binary_velocities("SHELL", shell_part, force_env%root_section, para_env, &
752 1794 : subsys_section, shellvel_read, cell)
753 : CALL read_binary_velocities("CORE", core_part, force_env%root_section, para_env, &
754 1794 : subsys_section, corevel_read, cell)
755 :
756 1794 : 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 1794 : 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 1528 : SUBROUTINE generate_velocities(simpar, part, force_env, globenv, md_env, &
1176 1528 : 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 1528 : NULLIFY (atomic_kind)
1192 1528 : natoms = SIZE(part)
1193 :
1194 437432 : DO i = 1, natoms
1195 435904 : atomic_kind => part(i)%atomic_kind
1196 435904 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1197 435904 : part(i)%v(1) = 0.0_dp
1198 435904 : part(i)%v(2) = 0.0_dp
1199 435904 : part(i)%v(3) = 0.0_dp
1200 437432 : IF (mass /= 0.0) THEN
1201 435970 : 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 430292 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1219 430292 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1220 865750 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1221 : END SELECT
1222 : END IF
1223 : END DO
1224 :
1225 1528 : CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1226 1528 : 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 1528 : IF (shell_present) THEN
1231 84 : CALL optimize_shell_core(force_env, part, shell_part, core_part, globenv)
1232 : END IF
1233 1528 : 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 1528 : 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 1528 : IF (simpar%soften_nsteps <= 0) RETURN !nothing todo
1259 :
1260 308 : IF (ANY(is_fixed /= use_perd_none)) THEN
1261 0 : CPABORT("Velocitiy softening with constraints is not supported.")
1262 : END IF
1263 :
1264 : !backup positions
1265 308 : DO i = 1, SIZE(part)
1266 1148 : x0(i, :) = part(i)%r
1267 : END DO
1268 :
1269 588 : DO k = 1, simpar%soften_nsteps
1270 :
1271 : !use normalized velocities as displace direction
1272 6160 : DO i = 1, SIZE(part)
1273 22960 : N(i, :) = part(i)%v
1274 : END DO
1275 37520 : N = N/SQRT(SUM(N**2))
1276 :
1277 : ! displace system temporarly to calculate forces
1278 6160 : DO i = 1, SIZE(part)
1279 22960 : part(i)%r = part(i)%r + simpar%soften_delta*N(i, :)
1280 : END DO
1281 560 : CALL force_env_calc_energy_force(force_env)
1282 :
1283 : ! calculate velocity update direction F_t
1284 6160 : DO i = 1, SIZE(part)
1285 22960 : F(i, :) = part(i)%f
1286 : END DO
1287 37520 : F_t = F - N*SUM(N*F)
1288 :
1289 : ! restore positions and update velocities
1290 6160 : DO i = 1, SIZE(part)
1291 22400 : part(i)%r = x0(i, :)
1292 22960 : part(i)%v = part(i)%v + simpar%soften_alpha*F_t(i, :)
1293 : END DO
1294 :
1295 588 : CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1296 : END DO
1297 :
1298 28 : IF (iw > 0) THEN
1299 0 : WRITE (iw, "(A,T71, I10)") " Velocities softening Steps: ", simpar%soften_nsteps
1300 0 : WRITE (iw, "(A,T71, E10.3)") " Velocities softening NORM(F_t): ", SQRT(SUM(F_t**2))
1301 : END IF
1302 28 : END SUBROUTINE soften_velocities
1303 :
1304 : ! **************************************************************************************************
1305 : !> \brief Scale velocities according to temperature and remove rigid body motion.
1306 : !> \param simpar ...
1307 : !> \param part ...
1308 : !> \param force_env ...
1309 : !> \param md_env ...
1310 : !> \param is_fixed ...
1311 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1312 : ! **************************************************************************************************
1313 2088 : SUBROUTINE normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1314 : TYPE(simpar_type), POINTER :: simpar
1315 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1316 : TYPE(force_env_type), POINTER :: force_env
1317 : TYPE(md_environment_type), POINTER :: md_env
1318 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1319 :
1320 : REAL(KIND=dp) :: ekin
1321 : REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
1322 : TYPE(cell_type), POINTER :: cell
1323 :
1324 2088 : NULLIFY (cell)
1325 :
1326 : ! Subtract the vcom
1327 2088 : CALL compute_vcom(part, is_fixed, vcom)
1328 2088 : CALL subtract_vcom(part, is_fixed, vcom)
1329 : ! If requested and the system is not periodic, subtract the angular velocity
1330 2088 : CALL force_env_get(force_env, cell=cell)
1331 8352 : IF (SUM(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN
1332 4 : CALL compute_rcom(part, is_fixed, rcom)
1333 4 : CALL compute_vang(part, is_fixed, rcom, vang)
1334 4 : CALL subtract_vang(part, is_fixed, rcom, vang)
1335 : END IF
1336 : ! Rescale the velocities
1337 2088 : IF (simpar%do_thermal_region) THEN
1338 2 : CALL rescale_vel_region(part, md_env, simpar)
1339 : ELSE
1340 2086 : ekin = compute_ekin(part)
1341 2086 : CALL rescale_vel(part, simpar, ekin)
1342 : END IF
1343 2088 : END SUBROUTINE normalize_velocities
1344 :
1345 : ! **************************************************************************************************
1346 : !> \brief Computes Ekin, VCOM and Temp for particles
1347 : !> \param subsys ...
1348 : !> \param md_ener ...
1349 : !> \param vsubtract ...
1350 : !> \par History
1351 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1352 : ! **************************************************************************************************
1353 42 : SUBROUTINE reset_vcom(subsys, md_ener, vsubtract)
1354 : TYPE(cp_subsys_type), POINTER :: subsys
1355 : TYPE(md_ener_type), POINTER :: md_ener
1356 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vsubtract
1357 :
1358 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_vcom'
1359 :
1360 : INTEGER :: atom, handle, iatom, ikind, natom, &
1361 : shell_index
1362 42 : INTEGER, DIMENSION(:), POINTER :: atom_list
1363 : LOGICAL :: is_shell
1364 : REAL(KIND=dp) :: ekin_old, imass_c, imass_s, mass, v2
1365 : REAL(KIND=dp), DIMENSION(3) :: tmp, v
1366 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1367 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1368 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1369 : shell_particles
1370 : TYPE(shell_kind_type), POINTER :: shell
1371 :
1372 42 : NULLIFY (particles, atomic_kind, atomic_kinds, atom_list, shell)
1373 42 : CALL timeset(routineN, handle)
1374 :
1375 : CALL cp_subsys_get(subsys, &
1376 : atomic_kinds=atomic_kinds, &
1377 : particles=particles, &
1378 : shell_particles=shell_particles, &
1379 42 : core_particles=core_particles)
1380 :
1381 42 : ekin_old = md_ener%ekin
1382 : ! Possibly subtract a quantity from all velocities
1383 126 : DO ikind = 1, atomic_kinds%n_els
1384 84 : atomic_kind => atomic_kinds%els(ikind)
1385 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, &
1386 84 : natom=natom, mass=mass, shell_active=is_shell, shell=shell)
1387 126 : IF (is_shell) THEN
1388 336 : tmp = 0.5_dp*vsubtract*mass
1389 84 : imass_s = 1.0_dp/shell%mass_shell
1390 84 : imass_c = 1.0_dp/shell%mass_core
1391 3780 : DO iatom = 1, natom
1392 3696 : atom = atom_list(iatom)
1393 3696 : shell_index = particles%els(atom)%shell_index
1394 14784 : shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s
1395 14784 : core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c
1396 14868 : particles%els(atom)%v = particles%els(atom)%v - vsubtract
1397 : END DO
1398 : ELSE
1399 0 : DO iatom = 1, natom
1400 0 : atom = atom_list(iatom)
1401 0 : particles%els(atom)%v = particles%els(atom)%v - vsubtract
1402 : END DO
1403 : END IF
1404 : END DO
1405 : ! Compute Kinetic Energy and COM Velocity
1406 168 : md_ener%vcom = 0.0_dp
1407 42 : md_ener%total_mass = 0.0_dp
1408 42 : md_ener%ekin = 0.0_dp
1409 126 : DO ikind = 1, atomic_kinds%n_els
1410 84 : atomic_kind => atomic_kinds%els(ikind)
1411 84 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom)
1412 84 : v2 = 0.0_dp
1413 84 : v = 0.0_dp
1414 3780 : DO iatom = 1, natom
1415 3696 : atom = atom_list(iatom)
1416 14784 : v2 = v2 + SUM(particles%els(atom)%v**2)
1417 3696 : v(1) = v(1) + particles%els(atom)%v(1)
1418 3696 : v(2) = v(2) + particles%els(atom)%v(2)
1419 3780 : v(3) = v(3) + particles%els(atom)%v(3)
1420 : END DO
1421 84 : md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2
1422 84 : md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1)
1423 84 : md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2)
1424 84 : md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3)
1425 210 : md_ener%total_mass = md_ener%total_mass + REAL(natom, KIND=dp)*mass
1426 : END DO
1427 168 : md_ener%vcom = md_ener%vcom/md_ener%total_mass
1428 42 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1429 42 : IF (md_ener%nfree /= 0) THEN
1430 42 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1431 : END IF
1432 42 : CALL timestop(handle)
1433 :
1434 42 : END SUBROUTINE reset_vcom
1435 :
1436 : ! **************************************************************************************************
1437 : !> \brief Scale velocities to get the correct temperature
1438 : !> \param subsys ...
1439 : !> \param md_ener ...
1440 : !> \param temp_expected ...
1441 : !> \param temp_tol ...
1442 : !> \param iw ...
1443 : !> \par History
1444 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1445 : ! **************************************************************************************************
1446 12914 : SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
1447 : TYPE(cp_subsys_type), POINTER :: subsys
1448 : TYPE(md_ener_type), POINTER :: md_ener
1449 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1450 : INTEGER, INTENT(IN) :: iw
1451 :
1452 : REAL(KIND=dp) :: ekin_old, scale, temp_old
1453 :
1454 12914 : IF (ABS(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN
1455 2640 : scale = 0.0_dp
1456 2640 : IF (md_ener%temp_part > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_part)*kelvin)
1457 2640 : ekin_old = md_ener%ekin
1458 2640 : temp_old = md_ener%temp_part
1459 2640 : md_ener%ekin = 0.0_dp
1460 2640 : md_ener%temp_part = 0.0_dp
1461 10560 : md_ener%vcom = 0.0_dp
1462 2640 : md_ener%total_mass = 0.0_dp
1463 :
1464 2640 : CALL scale_velocity_low(subsys, scale, ireg=0, ekin=md_ener%ekin, vcom=md_ener%vcom)
1465 2640 : IF (md_ener%nfree /= 0) THEN
1466 2640 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1467 : END IF
1468 2640 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1469 2640 : IF (iw > 0) THEN
1470 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1471 1320 : 'MD_VEL| Temperature scaled to requested temperature'
1472 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1473 1320 : 'MD_VEL| Old temperature [K]', temp_old, &
1474 2640 : 'MD_VEL| New temperature [K]', md_ener%temp_part
1475 : END IF
1476 : END IF
1477 :
1478 12914 : END SUBROUTINE scale_velocity
1479 :
1480 : ! **************************************************************************************************
1481 : !> \brief Scale velocities of set of regions
1482 : !> \param md_env ...
1483 : !> \param subsys ...
1484 : !> \param md_ener ...
1485 : !> \param simpar ...
1486 : !> \param iw ...
1487 : !> \par author MI
1488 : ! **************************************************************************************************
1489 96 : SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1490 :
1491 : TYPE(md_environment_type), POINTER :: md_env
1492 : TYPE(cp_subsys_type), POINTER :: subsys
1493 : TYPE(md_ener_type), POINTER :: md_ener
1494 : TYPE(simpar_type), POINTER :: simpar
1495 : INTEGER, INTENT(IN) :: iw
1496 :
1497 : INTEGER :: ireg, nfree, nfree_done, nregions
1498 : REAL(KIND=dp) :: ekin, ekin_old, ekin_total_new, fscale, &
1499 : vcom(3), vcom_total(3)
1500 96 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp_new, temp_old
1501 : TYPE(particle_list_type), POINTER :: particles
1502 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1503 : TYPE(thermal_region_type), POINTER :: t_region
1504 : TYPE(thermal_regions_type), POINTER :: thermal_regions
1505 :
1506 96 : NULLIFY (particles, part, thermal_regions, t_region)
1507 96 : CALL cp_subsys_get(subsys, particles=particles)
1508 96 : part => particles%els
1509 96 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
1510 :
1511 96 : nregions = thermal_regions%nregions
1512 96 : nfree_done = 0
1513 96 : ekin_total_new = 0.0_dp
1514 96 : ekin_old = md_ener%ekin
1515 : vcom_total = 0.0_dp
1516 384 : ALLOCATE (temp_new(0:nregions), temp_old(0:nregions))
1517 96 : temp_new = 0.0_dp
1518 96 : temp_old = 0.0_dp
1519 : !loop regions
1520 256 : DO ireg = 1, nregions
1521 160 : NULLIFY (t_region)
1522 160 : t_region => thermal_regions%thermal_region(ireg)
1523 160 : nfree = 3*t_region%npart
1524 160 : ekin = compute_ekin(part, ireg)
1525 160 : IF (nfree > 0) t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1526 160 : temp_old(ireg) = t_region%temperature
1527 160 : IF (t_region%temp_tol > 0.0_dp .AND. &
1528 : ABS(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN
1529 2 : fscale = SQRT((t_region%temp_expected/t_region%temperature)*kelvin)
1530 2 : CALL scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1531 2 : t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1532 2 : temp_new(ireg) = t_region%temperature
1533 : END IF
1534 160 : nfree_done = nfree_done + nfree
1535 256 : ekin_total_new = ekin_total_new + ekin
1536 : END DO
1537 96 : nfree = simpar%nfree - nfree_done
1538 96 : ekin = compute_ekin(part, ireg=0)
1539 96 : IF (nfree > 0) thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1540 96 : temp_old(0) = thermal_regions%temp_reg0
1541 96 : IF (simpar%temp_tol > 0.0_dp .AND. nfree > 0) THEN
1542 0 : IF (ABS(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN
1543 0 : fscale = SQRT((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin)
1544 0 : CALL scale_velocity_low(subsys, fscale, 0, ekin, vcom)
1545 0 : thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1546 0 : temp_new(0) = thermal_regions%temp_reg0
1547 : END IF
1548 : END IF
1549 96 : ekin_total_new = ekin_total_new + ekin
1550 :
1551 96 : md_ener%ekin = ekin_total_new
1552 96 : IF (md_ener%nfree /= 0) THEN
1553 96 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1554 : END IF
1555 96 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1556 96 : IF (iw > 0) THEN
1557 176 : DO ireg = 0, nregions
1558 176 : IF (temp_new(ireg) > 0.0_dp) THEN
1559 : WRITE (UNIT=iw, FMT='(/,T2,A,I0,A)') &
1560 1 : 'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
1561 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1562 1 : 'MD_VEL| Old temperature [K]', temp_old(ireg), &
1563 2 : 'MD_VEL| New temperature [K]', temp_new(ireg)
1564 : END IF
1565 : END DO
1566 : END IF
1567 96 : DEALLOCATE (temp_new, temp_old)
1568 :
1569 96 : END SUBROUTINE scale_velocity_region
1570 :
1571 : ! **************************************************************************************************
1572 : !> \brief Scale velocities for a specific region
1573 : !> \param subsys ...
1574 : !> \param fscale ...
1575 : !> \param ireg ...
1576 : !> \param ekin ...
1577 : !> \param vcom ...
1578 : !> \par author MI
1579 : ! **************************************************************************************************
1580 2642 : SUBROUTINE scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1581 :
1582 : TYPE(cp_subsys_type), POINTER :: subsys
1583 : REAL(KIND=dp), INTENT(IN) :: fscale
1584 : INTEGER, INTENT(IN) :: ireg
1585 : REAL(KIND=dp), INTENT(OUT) :: ekin, vcom(3)
1586 :
1587 : INTEGER :: atom, iatom, ikind, my_ireg, natom, &
1588 : shell_index
1589 2642 : INTEGER, DIMENSION(:), POINTER :: atom_list
1590 : LOGICAL :: is_shell
1591 : REAL(KIND=dp) :: imass, mass, tmass, v2
1592 : REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1593 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1594 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1595 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1596 : shell_particles
1597 : TYPE(shell_kind_type), POINTER :: shell
1598 :
1599 2642 : NULLIFY (atomic_kinds, particles, shell_particles, core_particles, shell, atom_list)
1600 :
1601 2642 : my_ireg = ireg
1602 2642 : ekin = 0.0_dp
1603 2642 : tmass = 0.0_dp
1604 2642 : vcom = 0.0_dp
1605 :
1606 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
1607 2642 : shell_particles=shell_particles, core_particles=core_particles)
1608 :
1609 9082 : DO ikind = 1, atomic_kinds%n_els
1610 6440 : atomic_kind => atomic_kinds%els(ikind)
1611 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, &
1612 6440 : natom=natom, shell_active=is_shell, shell=shell)
1613 6440 : IF (is_shell) THEN
1614 124 : imass = 1.0_dp/mass
1615 124 : v2 = 0.0_dp
1616 124 : v = 0.0_dp
1617 5740 : DO iatom = 1, natom
1618 5616 : atom = atom_list(iatom)
1619 : !check region
1620 5616 : IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
1621 :
1622 21080 : particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1623 5270 : shell_index = particles%els(atom)%shell_index
1624 21080 : vs = shell_particles%els(shell_index)%v
1625 21080 : vc = core_particles%els(shell_index)%v
1626 5270 : tmp(1) = imass*(vs(1) - vc(1))
1627 5270 : tmp(2) = imass*(vs(2) - vc(2))
1628 5270 : tmp(3) = imass*(vs(3) - vc(3))
1629 :
1630 5270 : shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core
1631 5270 : shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core
1632 5270 : shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core
1633 :
1634 5270 : core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell
1635 5270 : core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell
1636 5270 : core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell
1637 :
1638 : ! kinetic energy and velocity of COM
1639 21080 : v2 = v2 + SUM(particles%els(atom)%v**2)
1640 5270 : v(1) = v(1) + particles%els(atom)%v(1)
1641 5270 : v(2) = v(2) + particles%els(atom)%v(2)
1642 5270 : v(3) = v(3) + particles%els(atom)%v(3)
1643 5740 : tmass = tmass + mass
1644 : END DO
1645 : ELSE
1646 6316 : v2 = 0.0_dp
1647 6316 : v = 0.0_dp
1648 22826 : DO iatom = 1, natom
1649 16510 : atom = atom_list(iatom)
1650 : !check region
1651 16510 : IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
1652 :
1653 66040 : particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1654 : ! kinetic energy and velocity of COM
1655 66040 : v2 = v2 + SUM(particles%els(atom)%v**2)
1656 16510 : v(1) = v(1) + particles%els(atom)%v(1)
1657 16510 : v(2) = v(2) + particles%els(atom)%v(2)
1658 16510 : v(3) = v(3) + particles%els(atom)%v(3)
1659 22826 : tmass = tmass + mass
1660 : END DO
1661 : END IF
1662 6440 : ekin = ekin + 0.5_dp*mass*v2
1663 6440 : vcom(1) = vcom(1) + mass*v(1)
1664 6440 : vcom(2) = vcom(2) + mass*v(2)
1665 15522 : vcom(3) = vcom(3) + mass*v(3)
1666 :
1667 : END DO
1668 10568 : vcom = vcom/tmass
1669 :
1670 2642 : END SUBROUTINE scale_velocity_low
1671 :
1672 : ! **************************************************************************************************
1673 : !> \brief Scale internal motion of CORE-SHELL model to the correct temperature
1674 : !> \param subsys ...
1675 : !> \param md_ener ...
1676 : !> \param temp_expected ...
1677 : !> \param temp_tol ...
1678 : !> \param iw ...
1679 : !> \par History
1680 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1681 : ! **************************************************************************************************
1682 1060 : SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
1683 : TYPE(cp_subsys_type), POINTER :: subsys
1684 : TYPE(md_ener_type), POINTER :: md_ener
1685 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1686 : INTEGER, INTENT(IN) :: iw
1687 :
1688 : INTEGER :: atom, iatom, ikind, natom, shell_index
1689 1060 : INTEGER, DIMENSION(:), POINTER :: atom_list
1690 : LOGICAL :: is_shell
1691 : REAL(KIND=dp) :: ekin_shell_old, fac_mass, mass, scale, &
1692 : temp_shell_old, v2
1693 : REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1694 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1695 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1696 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1697 : shell_particles
1698 : TYPE(shell_kind_type), POINTER :: shell
1699 :
1700 1060 : NULLIFY (atom_list, atomic_kinds, atomic_kind, core_particles, particles, shell_particles, shell)
1701 1060 : IF (ABS(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN
1702 80 : scale = 0.0_dp
1703 80 : IF (md_ener%temp_shell > EPSILON(0.0_dp)) scale = SQRT((temp_expected/md_ener%temp_shell)*kelvin)
1704 80 : ekin_shell_old = md_ener%ekin_shell
1705 80 : temp_shell_old = md_ener%temp_shell
1706 80 : md_ener%ekin_shell = 0.0_dp
1707 80 : md_ener%temp_shell = 0.0_dp
1708 :
1709 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles, &
1710 80 : core_particles=core_particles)
1711 :
1712 240 : DO ikind = 1, atomic_kinds%n_els
1713 160 : atomic_kind => atomic_kinds%els(ikind)
1714 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom, &
1715 160 : shell_active=is_shell, shell=shell)
1716 240 : IF (is_shell) THEN
1717 160 : fac_mass = 1.0_dp/mass
1718 160 : v2 = 0.0_dp
1719 776 : DO iatom = 1, natom
1720 616 : atom = atom_list(iatom)
1721 616 : shell_index = particles%els(atom)%shell_index
1722 2464 : vs = shell_particles%els(shell_index)%v
1723 2464 : vc = core_particles%els(shell_index)%v
1724 2464 : v = particles%els(atom)%v
1725 616 : tmp(1) = fac_mass*(vc(1) - vs(1))
1726 616 : tmp(2) = fac_mass*(vc(2) - vs(2))
1727 616 : tmp(3) = fac_mass*(vc(3) - vs(3))
1728 :
1729 616 : shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1)
1730 616 : shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2)
1731 616 : shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3)
1732 :
1733 616 : core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1)
1734 616 : core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2)
1735 616 : core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3)
1736 :
1737 2464 : vs = shell_particles%els(shell_index)%v
1738 2464 : vc = core_particles%els(shell_index)%v
1739 616 : tmp(1) = vc(1) - vs(1)
1740 616 : tmp(2) = vc(2) - vs(2)
1741 616 : tmp(3) = vc(3) - vs(3)
1742 2624 : v2 = v2 + SUM(tmp**2)
1743 : END DO
1744 160 : md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2
1745 : END IF
1746 : END DO
1747 80 : IF (md_ener%nfree_shell > 0) THEN
1748 80 : md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)*kelvin
1749 : END IF
1750 80 : md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell
1751 80 : IF (iw > 0) THEN
1752 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1753 40 : 'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
1754 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1755 40 : 'MD_VEL| Old temperature [K]', temp_shell_old, &
1756 80 : 'MD_VEL| New temperature [K]', md_ener%temp_shell
1757 : END IF
1758 : END IF
1759 :
1760 1060 : END SUBROUTINE scale_velocity_internal
1761 :
1762 : ! **************************************************************************************************
1763 : !> \brief Scale barostat velocities to get the desired temperature
1764 : !> \param md_env ...
1765 : !> \param md_ener ...
1766 : !> \param temp_expected ...
1767 : !> \param temp_tol ...
1768 : !> \param iw ...
1769 : !> \par History
1770 : !> MI 02.2008
1771 : ! **************************************************************************************************
1772 40 : SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
1773 : TYPE(md_environment_type), POINTER :: md_env
1774 : TYPE(md_ener_type), POINTER :: md_ener
1775 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1776 : INTEGER, INTENT(IN) :: iw
1777 :
1778 : INTEGER :: i, j, nfree
1779 : REAL(KIND=dp) :: ekin_old, scale, temp_old
1780 40 : TYPE(npt_info_type), POINTER :: npt(:, :)
1781 : TYPE(simpar_type), POINTER :: simpar
1782 :
1783 40 : NULLIFY (npt, simpar)
1784 40 : CALL get_md_env(md_env, simpar=simpar, npt=npt)
1785 40 : IF (ABS(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN
1786 2 : scale = 0.0_dp
1787 2 : IF (md_ener%temp_baro > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_baro)*kelvin)
1788 2 : ekin_old = md_ener%baro_kin
1789 2 : temp_old = md_ener%temp_baro
1790 2 : md_ener%baro_kin = 0.0_dp
1791 2 : md_ener%temp_baro = 0.0_dp
1792 : IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
1793 2 : .OR. simpar%ensemble == npt_ia_ensemble) THEN
1794 0 : npt(1, 1)%v = npt(1, 1)%v*scale
1795 0 : md_ener%baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
1796 : ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
1797 2 : md_ener%baro_kin = 0.0_dp
1798 8 : DO i = 1, 3
1799 26 : DO j = 1, 3
1800 18 : npt(i, j)%v = npt(i, j)%v*scale
1801 24 : md_ener%baro_kin = md_ener%baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
1802 : END DO
1803 : END DO
1804 : END IF
1805 2 : nfree = SIZE(npt, 1)*SIZE(npt, 2)
1806 2 : md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/REAL(nfree, dp)*kelvin
1807 2 : IF (iw > 0) THEN
1808 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1809 1 : 'MD_VEL| Temperature of barostat motion scaled to requested temperature'
1810 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1811 1 : 'MD_VEL| Old temperature [K]', temp_old, &
1812 2 : 'MD_VEL| New temperature [K]', md_ener%temp_baro
1813 : END IF
1814 : END IF
1815 :
1816 40 : END SUBROUTINE scale_velocity_baro
1817 :
1818 : ! **************************************************************************************************
1819 : !> \brief Perform all temperature manipulations during a QS MD run.
1820 : !> \param simpar ...
1821 : !> \param md_env ...
1822 : !> \param md_ener ...
1823 : !> \param force_env ...
1824 : !> \param logger ...
1825 : !> \par History
1826 : !> Creation (15.09.2003,MK)
1827 : !> adapted to force_env (05.10.2003,fawzi)
1828 : !> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1829 : ! **************************************************************************************************
1830 40643 : SUBROUTINE temperature_control(simpar, md_env, md_ener, force_env, logger)
1831 :
1832 : TYPE(simpar_type), POINTER :: simpar
1833 : TYPE(md_environment_type), POINTER :: md_env
1834 : TYPE(md_ener_type), POINTER :: md_ener
1835 : TYPE(force_env_type), POINTER :: force_env
1836 : TYPE(cp_logger_type), POINTER :: logger
1837 :
1838 : CHARACTER(LEN=*), PARAMETER :: routineN = 'temperature_control'
1839 :
1840 : INTEGER :: handle, iw
1841 : TYPE(cp_subsys_type), POINTER :: subsys
1842 : TYPE(mp_para_env_type), POINTER :: para_env
1843 :
1844 40643 : CALL timeset(routineN, handle)
1845 40643 : NULLIFY (subsys, para_env)
1846 40643 : CPASSERT(ASSOCIATED(simpar))
1847 40643 : CPASSERT(ASSOCIATED(md_ener))
1848 40643 : CPASSERT(ASSOCIATED(force_env))
1849 40643 : CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
1850 : iw = cp_print_key_unit_nr(logger, force_env%root_section, "MOTION%MD%PRINT%PROGRAM_RUN_INFO", &
1851 40643 : extension=".mdLog")
1852 :
1853 : ! Control the particle motion
1854 40643 : IF (simpar%do_thermal_region) THEN
1855 96 : CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1856 : ELSE
1857 40547 : IF (simpar%temp_tol > 0.0_dp) THEN
1858 12870 : CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw)
1859 : END IF
1860 : END IF
1861 : ! Control the internal core-shell motion
1862 40643 : IF (simpar%temp_sh_tol > 0.0_dp) THEN
1863 1060 : CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw)
1864 : END IF
1865 : ! Control cell motion
1866 43183 : SELECT CASE (simpar%ensemble)
1867 : CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
1868 : npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_ia_ensemble)
1869 40643 : IF (simpar%temp_baro_tol > 0.0_dp) THEN
1870 40 : CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw)
1871 : END IF
1872 : END SELECT
1873 :
1874 : CALL cp_print_key_finished_output(iw, logger, force_env%root_section, &
1875 40643 : "MOTION%MD%PRINT%PROGRAM_RUN_INFO")
1876 40643 : CALL timestop(handle)
1877 40643 : END SUBROUTINE temperature_control
1878 :
1879 : ! **************************************************************************************************
1880 : !> \brief Set to 0 the velocity of the COM along MD runs, if required.
1881 : !> \param md_ener ...
1882 : !> \param force_env ...
1883 : !> \param md_section ...
1884 : !> \param logger ...
1885 : !> \par History
1886 : !> Creation (29.04.2007,MI)
1887 : !> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1888 : ! **************************************************************************************************
1889 81286 : SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
1890 :
1891 : TYPE(md_ener_type), POINTER :: md_ener
1892 : TYPE(force_env_type), POINTER :: force_env
1893 : TYPE(section_vals_type), POINTER :: md_section
1894 : TYPE(cp_logger_type), POINTER :: logger
1895 :
1896 : CHARACTER(LEN=*), PARAMETER :: routineN = 'comvel_control'
1897 :
1898 : INTEGER :: handle, iw
1899 : LOGICAL :: explicit
1900 : REAL(KIND=dp) :: comvel_tol, temp_old, vel_com
1901 : REAL(KIND=dp), DIMENSION(3) :: vcom_old
1902 : TYPE(cp_subsys_type), POINTER :: subsys
1903 :
1904 40643 : CALL timeset(routineN, handle)
1905 40643 : NULLIFY (subsys)
1906 40643 : CPASSERT(ASSOCIATED(force_env))
1907 40643 : CALL force_env_get(force_env, subsys=subsys)
1908 :
1909 : ! Print COMVEL and COM Position
1910 40643 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
1911 40643 : IF (iw > 0) THEN
1912 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1913 7785 : "MD_VEL| Centre of mass motion (COM)"
1914 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1915 31140 : "MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
1916 : END IF
1917 40643 : CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")
1918 :
1919 : ! If requested rescale COMVEL
1920 40643 : CALL section_vals_val_get(md_section, "COMVEL_TOL", explicit=explicit)
1921 40643 : IF (explicit) THEN
1922 826 : CALL section_vals_val_get(md_section, "COMVEL_TOL", r_val=comvel_tol)
1923 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1924 826 : extension=".mdLog")
1925 826 : vel_com = SQRT(md_ener%vcom(1)**2 + md_ener%vcom(2)**2 + md_ener%vcom(3)**2)
1926 :
1927 : ! Subtract the velocity of the COM, if requested
1928 826 : IF (vel_com > comvel_tol) THEN
1929 42 : temp_old = md_ener%temp_part/kelvin
1930 168 : vcom_old = md_ener%vcom
1931 42 : CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
1932 42 : CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
1933 42 : IF (iw > 0) THEN
1934 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1935 21 : "MD_VEL| Old VCOM [a.u.]", vcom_old(1:3)
1936 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1937 84 : "MD_VEL| New VCOM [a.u.]", md_ener%vcom(1:3)
1938 : END IF
1939 : END IF
1940 : CALL cp_print_key_finished_output(iw, logger, md_section, &
1941 826 : "PRINT%PROGRAM_RUN_INFO")
1942 : END IF
1943 :
1944 40643 : CALL timestop(handle)
1945 40643 : END SUBROUTINE comvel_control
1946 :
1947 : ! **************************************************************************************************
1948 : !> \brief Set to 0 the angular velocity along MD runs, if required.
1949 : !> \param md_ener ...
1950 : !> \param force_env ...
1951 : !> \param md_section ...
1952 : !> \param logger ...
1953 : !> \par History
1954 : !> Creation (10.2009) Teodoro Laino [tlaino]
1955 : ! **************************************************************************************************
1956 40643 : SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
1957 :
1958 : TYPE(md_ener_type), POINTER :: md_ener
1959 : TYPE(force_env_type), POINTER :: force_env
1960 : TYPE(section_vals_type), POINTER :: md_section
1961 : TYPE(cp_logger_type), POINTER :: logger
1962 :
1963 : CHARACTER(LEN=*), PARAMETER :: routineN = 'angvel_control'
1964 :
1965 : INTEGER :: handle, ifixd, imolecule_kind, iw, natoms
1966 40643 : INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
1967 : LOGICAL :: explicit
1968 : REAL(KIND=dp) :: angvel_tol, rcom(3), temp_old, vang(3), &
1969 : vang_new(3)
1970 : TYPE(cell_type), POINTER :: cell
1971 : TYPE(cp_subsys_type), POINTER :: subsys
1972 40643 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
1973 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1974 40643 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1975 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1976 : TYPE(particle_list_type), POINTER :: particles
1977 :
1978 40643 : CALL timeset(routineN, handle)
1979 : ! If requested rescale ANGVEL
1980 40643 : CALL section_vals_val_get(md_section, "ANGVEL_TOL", explicit=explicit)
1981 40643 : IF (explicit) THEN
1982 40 : NULLIFY (subsys, cell)
1983 40 : CPASSERT(ASSOCIATED(force_env))
1984 40 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
1985 :
1986 160 : IF (SUM(cell%perd(1:3)) == 0) THEN
1987 40 : CALL section_vals_val_get(md_section, "ANGVEL_TOL", r_val=angvel_tol)
1988 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1989 40 : extension=".mdLog")
1990 :
1991 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, &
1992 40 : particles=particles)
1993 :
1994 40 : natoms = SIZE(particles%els)
1995 : ! Build a list of all fixed atoms (if any)
1996 120 : ALLOCATE (is_fixed(natoms))
1997 :
1998 40 : is_fixed = use_perd_none
1999 40 : molecule_kind_set => molecule_kinds%els
2000 600 : DO imolecule_kind = 1, molecule_kinds%n_els
2001 560 : molecule_kind => molecule_kind_set(imolecule_kind)
2002 560 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2003 600 : IF (ASSOCIATED(fixd_list)) THEN
2004 0 : DO ifixd = 1, SIZE(fixd_list)
2005 0 : IF (.NOT. fixd_list(ifixd)%restraint%active) THEN
2006 0 : is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2007 : END IF
2008 : END DO
2009 : END IF
2010 : END DO
2011 :
2012 : ! If requested and the system is not periodic, subtract the angular velocity
2013 40 : CALL compute_rcom(particles%els, is_fixed, rcom)
2014 40 : CALL compute_vang(particles%els, is_fixed, rcom, vang)
2015 : ! SQRT(DOT_PRODUCT(vang,vang))>angvel_tol
2016 160 : IF (DOT_PRODUCT(vang, vang) > (angvel_tol*angvel_tol)) THEN
2017 2 : CALL subtract_vang(particles%els, is_fixed, rcom, vang)
2018 :
2019 : ! Rescale velocities after removal
2020 2 : temp_old = md_ener%temp_part/kelvin
2021 2 : CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
2022 2 : CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
2023 2 : IF (iw > 0) THEN
2024 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
2025 1 : 'MD_VEL| Old VANG [a.u.]', vang(1:3)
2026 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
2027 1 : 'MD_VEL| New VANG [a.u.]', vang_new(1:3)
2028 : END IF
2029 : END IF
2030 :
2031 40 : DEALLOCATE (is_fixed)
2032 :
2033 : CALL cp_print_key_finished_output(iw, logger, md_section, &
2034 80 : "PRINT%PROGRAM_RUN_INFO")
2035 : END IF
2036 : END IF
2037 :
2038 40643 : CALL timestop(handle)
2039 81286 : END SUBROUTINE angvel_control
2040 :
2041 : ! **************************************************************************************************
2042 : !> \brief Initialize Velocities for MD runs
2043 : !> \param force_env ...
2044 : !> \param simpar ...
2045 : !> \param globenv ...
2046 : !> \param md_env ...
2047 : !> \param md_section ...
2048 : !> \param constraint_section ...
2049 : !> \param write_binary_restart_file ...
2050 : !> \par History
2051 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
2052 : ! **************************************************************************************************
2053 3540 : SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, &
2054 : constraint_section, write_binary_restart_file)
2055 :
2056 : TYPE(force_env_type), POINTER :: force_env
2057 : TYPE(simpar_type), POINTER :: simpar
2058 : TYPE(global_environment_type), POINTER :: globenv
2059 : TYPE(md_environment_type), POINTER :: md_env
2060 : TYPE(section_vals_type), POINTER :: md_section, constraint_section
2061 : LOGICAL, INTENT(IN) :: write_binary_restart_file
2062 :
2063 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_velocities'
2064 :
2065 : INTEGER :: handle, nconstraint, nconstraint_fixd
2066 : LOGICAL :: apply_cns0, shell_adiabatic, &
2067 : shell_present
2068 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2069 : TYPE(cell_type), POINTER :: cell
2070 : TYPE(cp_subsys_type), POINTER :: subsys
2071 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2072 : TYPE(mp_para_env_type), POINTER :: para_env
2073 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
2074 : shell_particles
2075 1770 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
2076 1770 : shell_particle_set
2077 : TYPE(section_vals_type), POINTER :: force_env_section, print_section, &
2078 : subsys_section
2079 :
2080 1770 : CALL timeset(routineN, handle)
2081 :
2082 1770 : NULLIFY (atomic_kinds, cell, para_env, subsys, molecule_kinds, core_particles, particles)
2083 1770 : NULLIFY (shell_particles, core_particle_set, particle_set, shell_particle_set)
2084 1770 : NULLIFY (force_env_section, print_section, subsys_section)
2085 :
2086 1770 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
2087 1770 : apply_cns0 = .FALSE.
2088 1770 : IF (simpar%constraint) THEN
2089 310 : CALL section_vals_val_get(constraint_section, "CONSTRAINT_INIT", l_val=apply_cns0)
2090 : END IF
2091 : ! Always initialize velocities and possibly restart them
2092 : CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env, &
2093 1770 : force_env_section=force_env_section)
2094 1770 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2095 :
2096 : CALL cp_subsys_get(subsys, &
2097 : atomic_kinds=atomic_kinds, &
2098 : core_particles=core_particles, &
2099 : molecule_kinds=molecule_kinds, &
2100 : particles=particles, &
2101 1770 : shell_particles=shell_particles)
2102 :
2103 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
2104 : shell_present=shell_present, &
2105 1770 : shell_adiabatic=shell_adiabatic)
2106 :
2107 1770 : NULLIFY (core_particle_set)
2108 : NULLIFY (particle_set)
2109 1770 : NULLIFY (shell_particle_set)
2110 1770 : particle_set => particles%els
2111 :
2112 1770 : IF (shell_present .AND. shell_adiabatic) THEN
2113 : ! Constraints are not yet implemented for core-shell models generally
2114 : CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els, &
2115 : nconstraint=nconstraint, &
2116 132 : nconstraint_fixd=nconstraint_fixd)
2117 132 : IF (nconstraint - nconstraint_fixd /= 0) THEN
2118 0 : CPABORT("Only the fixed atom constraint is implemented for core-shell models")
2119 : END IF
2120 : !MK CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,failure)
2121 132 : CPASSERT(ASSOCIATED(shell_particles))
2122 132 : CPASSERT(ASSOCIATED(core_particles))
2123 132 : shell_particle_set => shell_particles%els
2124 132 : core_particle_set => core_particles%els
2125 : END IF
2126 :
2127 : CALL initialize_velocities(simpar, &
2128 : particle_set, &
2129 : molecule_kinds=molecule_kinds, &
2130 : force_env=force_env, &
2131 : globenv=globenv, &
2132 : md_env=md_env, &
2133 : label="Velocities initialization", &
2134 : print_section=print_section, &
2135 : subsys_section=subsys_section, &
2136 : shell_present=(shell_present .AND. shell_adiabatic), &
2137 : shell_part=shell_particle_set, &
2138 : core_part=core_particle_set, &
2139 : force_rescaling=.FALSE., &
2140 : para_env=para_env, &
2141 3408 : write_binary_restart_file=write_binary_restart_file)
2142 :
2143 : ! Apply constraints if required and rescale velocities..
2144 1770 : IF (simpar%ensemble /= reftraj_ensemble) THEN
2145 1734 : IF (apply_cns0) THEN
2146 24 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
2147 : CALL force_env_shake(force_env, &
2148 : shake_tol=simpar%shake_tol, &
2149 : log_unit=simpar%info_constraint, &
2150 : lagrange_mult=simpar%lagrange_multipliers, &
2151 : dump_lm=simpar%dump_lm, &
2152 24 : compold=.TRUE.)
2153 : CALL force_env_rattle(force_env, shake_tol=simpar%shake_tol, &
2154 : log_unit=simpar%info_constraint, lagrange_mult=simpar%lagrange_multipliers, &
2155 24 : dump_lm=simpar%dump_lm, reset=.TRUE.)
2156 24 : IF (simpar%do_respa) THEN
2157 : CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
2158 0 : calc_force=.TRUE.)
2159 : CALL force_env_shake(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, compold=.TRUE.)
2162 : CALL force_env_rattle(force_env%sub_force_env(1)%force_env, &
2163 : shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2164 0 : lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, reset=.TRUE.)
2165 : END IF
2166 : ! Reinitialize velocities rescaling properly after rattle
2167 24 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2168 24 : CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
2169 : CALL initialize_velocities(simpar, &
2170 : particle_set, &
2171 : molecule_kinds=molecule_kinds, &
2172 : force_env=force_env, &
2173 : globenv=globenv, &
2174 : md_env=md_env, &
2175 : label="Re-Initializing velocities after applying constraints", &
2176 : print_section=print_section, &
2177 : subsys_section=subsys_section, &
2178 : shell_present=(shell_present .AND. shell_adiabatic), &
2179 : shell_part=shell_particle_set, &
2180 : core_part=core_particle_set, &
2181 : force_rescaling=.TRUE., &
2182 : para_env=para_env, &
2183 48 : write_binary_restart_file=write_binary_restart_file)
2184 : END IF
2185 : END IF
2186 :
2187 : ! Perform setup for a cascade run
2188 1770 : CALL initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2189 :
2190 1770 : CALL timestop(handle)
2191 :
2192 1770 : END SUBROUTINE setup_velocities
2193 :
2194 : ! **************************************************************************************************
2195 : !> \brief Perform the initialization for a cascade run
2196 : !> \param simpar ...
2197 : !> \param particle_set ...
2198 : !> \param molecule_kinds ...
2199 : !> \param md_section ...
2200 : !> \date 05.02.2012
2201 : !> \author Matthias Krack (MK)
2202 : !> \version 1.0
2203 : ! **************************************************************************************************
2204 1770 : SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2205 :
2206 : TYPE(simpar_type), POINTER :: simpar
2207 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2208 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2209 : TYPE(section_vals_type), POINTER :: md_section
2210 :
2211 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_cascade'
2212 :
2213 : CHARACTER(len=2*default_string_length) :: line
2214 : INTEGER :: handle, iatom, ifixd, imolecule_kind, &
2215 : iparticle, iw, natom, nparticle
2216 1770 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_index, is_fixed
2217 : LOGICAL :: init_cascade, is_ok, no_read_error
2218 : REAL(KIND=dp) :: ecom, ekin, energy, norm, temp, &
2219 : temperature
2220 1770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: matom, weight
2221 1770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vatom
2222 : REAL(KIND=dp), DIMENSION(3) :: vcom
2223 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2224 : TYPE(cp_logger_type), POINTER :: logger
2225 : TYPE(cp_sll_val_type), POINTER :: atom_list
2226 1770 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
2227 1770 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2228 : TYPE(molecule_kind_type), POINTER :: molecule_kind
2229 : TYPE(section_vals_type), POINTER :: atom_list_section, cascade_section, &
2230 : print_section
2231 : TYPE(val_type), POINTER :: val
2232 :
2233 1770 : CALL timeset(routineN, handle)
2234 :
2235 1770 : NULLIFY (atom_list)
2236 1770 : NULLIFY (atom_list_section)
2237 1770 : NULLIFY (atomic_kind)
2238 1770 : NULLIFY (cascade_section)
2239 1770 : NULLIFY (fixd_list)
2240 1770 : NULLIFY (molecule_kind)
2241 1770 : NULLIFY (molecule_kind_set)
2242 1770 : NULLIFY (logger)
2243 1770 : NULLIFY (val)
2244 :
2245 1770 : logger => cp_get_default_logger()
2246 1770 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
2247 1770 : iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
2248 :
2249 1770 : cascade_section => section_vals_get_subs_vals(md_section, "CASCADE")
2250 1770 : CALL section_vals_val_get(cascade_section, "_SECTION_PARAMETERS_", l_val=init_cascade)
2251 :
2252 1770 : nparticle = SIZE(particle_set)
2253 :
2254 1770 : IF (init_cascade) THEN
2255 :
2256 2 : CALL section_vals_val_get(cascade_section, "ENERGY", r_val=energy)
2257 2 : IF (energy < 0.0_dp) THEN
2258 0 : CPABORT("Error occurred reading &CASCADE section: Negative energy found")
2259 : END IF
2260 :
2261 2 : IF (iw > 0) THEN
2262 1 : ekin = cp_unit_from_cp2k(energy, "keV")
2263 : WRITE (UNIT=iw, FMT="(/,T2,A,T61,F20.6)") &
2264 1 : "CASCADE| Energy [keV]", ekin
2265 : WRITE (UNIT=iw, FMT="(T2,A)") &
2266 1 : "CASCADE|"
2267 : END IF
2268 :
2269 : ! Read the atomic velocities given in the input file
2270 2 : atom_list_section => section_vals_get_subs_vals(cascade_section, "ATOM_LIST")
2271 2 : CALL section_vals_val_get(atom_list_section, "_DEFAULT_KEYWORD_", n_rep_val=natom)
2272 2 : CALL section_vals_list_get(atom_list_section, "_DEFAULT_KEYWORD_", list=atom_list)
2273 2 : IF (natom <= 0) THEN
2274 0 : CPABORT("Error occurred reading &CASCADE section: No atom list found")
2275 : END IF
2276 :
2277 2 : IF (iw > 0) THEN
2278 : WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
2279 1 : "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "weight"
2280 : END IF
2281 :
2282 6 : ALLOCATE (atom_index(natom))
2283 6 : ALLOCATE (matom(natom))
2284 6 : ALLOCATE (vatom(3, natom))
2285 4 : ALLOCATE (weight(natom))
2286 :
2287 8 : DO iatom = 1, natom
2288 6 : is_ok = cp_sll_val_next(atom_list, val)
2289 6 : CALL val_get(val, c_val=line)
2290 : ! Read atomic index, velocity vector, and weight
2291 6 : no_read_error = .FALSE.
2292 6 : READ (UNIT=line, FMT=*, ERR=999) atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2293 : no_read_error = .TRUE.
2294 : 999 IF (.NOT. no_read_error) THEN
2295 0 : CPABORT("Error occurred reading &CASCADE section. Last line read <"//TRIM(line)//">")
2296 : END IF
2297 6 : IF ((atom_index(iatom) <= 0) .OR. ((atom_index(iatom) > nparticle))) THEN
2298 0 : CPABORT("Error occurred reading &CASCADE section: Invalid atom index found")
2299 : END IF
2300 6 : IF (weight(iatom) < 0.0_dp) THEN
2301 0 : CPABORT("Error occurred reading &CASCADE section: Negative weight found")
2302 : END IF
2303 8 : IF (iw > 0) THEN
2304 : WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
2305 3 : "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2306 : END IF
2307 : END DO
2308 :
2309 : ! Normalise velocities and weights
2310 : norm = 0.0_dp
2311 8 : DO iatom = 1, natom
2312 6 : iparticle = atom_index(iatom)
2313 6 : IF (particle_set(iparticle)%shell_index /= 0) THEN
2314 0 : CPWARN("Warning: The primary knock-on atom is a core-shell atom")
2315 : END IF
2316 6 : atomic_kind => particle_set(iparticle)%atomic_kind
2317 6 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
2318 8 : norm = norm + matom(iatom)*weight(iatom)
2319 : END DO
2320 8 : weight(:) = matom(:)*weight(:)*energy/norm
2321 8 : DO iatom = 1, natom
2322 24 : norm = NORM2(vatom(1:3, iatom))
2323 26 : vatom(1:3, iatom) = vatom(1:3, iatom)/norm
2324 : END DO
2325 :
2326 2 : IF (iw > 0) THEN
2327 : WRITE (UNIT=iw, FMT="(T2,A)") &
2328 1 : "CASCADE|", &
2329 1 : "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
2330 2 : "CASCADE|"
2331 : WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
2332 1 : "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
2333 4 : DO iatom = 1, natom
2334 3 : ekin = cp_unit_from_cp2k(weight(iatom), "keV")
2335 : WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
2336 4 : "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
2337 : END DO
2338 : END IF
2339 :
2340 : ! Apply velocity modifications
2341 8 : DO iatom = 1, natom
2342 6 : iparticle = atom_index(iatom)
2343 : particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
2344 26 : SQRT(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
2345 : END DO
2346 :
2347 2 : DEALLOCATE (atom_index)
2348 2 : DEALLOCATE (matom)
2349 2 : DEALLOCATE (vatom)
2350 2 : DEALLOCATE (weight)
2351 :
2352 6 : IF (iw > 0) THEN
2353 : ! Build a list of all fixed atoms (if any)
2354 3 : ALLOCATE (is_fixed(nparticle))
2355 1 : is_fixed = use_perd_none
2356 1 : molecule_kind_set => molecule_kinds%els
2357 2 : DO imolecule_kind = 1, molecule_kinds%n_els
2358 1 : molecule_kind => molecule_kind_set(imolecule_kind)
2359 1 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2360 2 : IF (ASSOCIATED(fixd_list)) THEN
2361 0 : DO ifixd = 1, SIZE(fixd_list)
2362 0 : IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2363 : END DO
2364 : END IF
2365 : END DO
2366 : ! Compute vcom, ecom and ekin for printout
2367 1 : CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
2368 1 : ekin = compute_ekin(particle_set) - ecom
2369 1 : IF (simpar%nfree == 0) THEN
2370 0 : CPASSERT(ekin == 0.0_dp)
2371 0 : temp = 0.0_dp
2372 : ELSE
2373 1 : temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
2374 : END IF
2375 1 : temperature = cp_unit_from_cp2k(temp, "K")
2376 : WRITE (UNIT=iw, FMT="(T2,A)") &
2377 1 : "CASCADE|"
2378 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.6)") &
2379 1 : "CASCADE| Temperature after cascade initialization [K]", temperature
2380 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,ES16.8))") &
2381 1 : "CASCADE| COM velocity", vcom(1:3)
2382 : !MK ! compute and log rcom and vang if not periodic
2383 : !MK CALL force_env_get(force_env,cell=cell)
2384 : !MK IF (SUM(cell%perd(1:3)) == 0) THEN
2385 : !MK CALL compute_rcom(particle_set,is_fixed,rcom)
2386 : !MK CALL compute_vang(particle_set,is_fixed,rcom,vang)
2387 : !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
2388 : !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
2389 : !MK END IF
2390 2 : DEALLOCATE (is_fixed)
2391 : END IF
2392 :
2393 : END IF
2394 :
2395 1770 : CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
2396 :
2397 1770 : CALL timestop(handle)
2398 :
2399 3540 : END SUBROUTINE initialize_cascade
2400 :
2401 : END MODULE md_vel_utils
|
