(git:97501a3)
Loading...
Searching...
No Matches
force_env_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Util force_env module
10!> \author Teodoro Laino [tlaino] - 02.2011
11! **************************************************************************************************
12MODULE force_env_utils
13
14 USE atomic_kind_list_types, ONLY: atomic_kind_list_type
15 USE cell_types, ONLY: cell_type
16 USE constraint, ONLY: rattle_control,&
17 shake_control
18 USE constraint_util, ONLY: getold
19 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
20 USE cp_subsys_types, ONLY: cp_subsys_get,&
21 cp_subsys_type
22 USE cp_units, ONLY: cp_unit_from_cp2k
23 USE distribution_1d_types, ONLY: distribution_1d_type
24 USE force_env_types, ONLY: force_env_get,&
25 force_env_type
26 USE input_section_types, ONLY: section_vals_get,&
27 section_vals_get_subs_vals,&
28 section_vals_type,&
29 section_vals_val_get
30 USE kinds, ONLY: default_string_length,&
31 dp
32 USE molecule_kind_list_types, ONLY: molecule_kind_list_type
33 USE molecule_list_types, ONLY: molecule_list_type
34 USE molecule_types, ONLY: global_constraint_type
35 USE particle_list_types, ONLY: particle_list_type
36 USE particle_types, ONLY: update_particle_set
37 USE physcon, ONLY: angstrom
38 USE string_utilities, ONLY: lowercase,&
39 uppercase
40#include "./base/base_uses.f90"
41
42 IMPLICIT NONE
43
44 PRIVATE
45
46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_utils'
47
48 PUBLIC :: force_env_shake, &
49 force_env_rattle, &
50 rescale_forces, &
51 write_atener, &
52 write_forces
53
54CONTAINS
55
56! **************************************************************************************************
57!> \brief perform shake (enforcing of constraints)
58!> \param force_env the force env to shake
59!> \param dt the dt for shake (if you are not interested in the velocities
60!> it can be any positive number)
61!> \param shake_tol the tolerance for shake
62!> \param log_unit if >0 then some information on the shake is printed,
63!> defaults to -1
64!> \param lagrange_mult ...
65!> \param dump_lm ...
66!> \param pos ...
67!> \param vel ...
68!> \param compold ...
69!> \param reset ...
70!> \author fawzi
71! **************************************************************************************************
72 SUBROUTINE force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
73 pos, vel, compold, reset)
74
75 TYPE(force_env_type), POINTER :: force_env
76 REAL(kind=dp), INTENT(IN), OPTIONAL :: dt
77 REAL(kind=dp), INTENT(IN) :: shake_tol
78 INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult
79 LOGICAL, INTENT(IN), OPTIONAL :: dump_lm
80 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
81 OPTIONAL, TARGET :: pos, vel
82 LOGICAL, INTENT(IN), OPTIONAL :: compold, reset
83
84 CHARACTER(len=*), PARAMETER :: routinen = 'force_env_shake'
85
86 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, &
87 my_log_unit, nparticle_kind, nparticle_local
88 LOGICAL :: has_pos, has_vel, my_dump_lm
89 REAL(kind=dp) :: mydt
90 REAL(kind=dp), DIMENSION(:, :), POINTER :: my_pos, my_vel
91 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
92 TYPE(cell_type), POINTER :: cell
93 TYPE(cp_subsys_type), POINTER :: subsys
94 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
95 TYPE(global_constraint_type), POINTER :: gci
96 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
97 TYPE(molecule_list_type), POINTER :: molecules
98 TYPE(particle_list_type), POINTER :: particles
99
100 CALL timeset(routinen, handle)
101 cpassert(ASSOCIATED(force_env))
102 cpassert(force_env%ref_count > 0)
103 my_log_unit = -1
104 IF (PRESENT(log_unit)) my_log_unit = log_unit
105 my_lagrange_mult = -1
106 IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
107 my_dump_lm = .false.
108 IF (PRESENT(dump_lm)) my_dump_lm = dump_lm
109 NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
110 my_pos, my_vel, gci)
111 IF (PRESENT(pos)) my_pos => pos
112 IF (PRESENT(vel)) my_vel => vel
113 mydt = 0.1_dp
114 IF (PRESENT(dt)) mydt = dt
115 CALL force_env_get(force_env, subsys=subsys, cell=cell)
116 CALL cp_subsys_get(subsys, &
117 atomic_kinds=atomic_kinds, &
118 local_molecules=local_molecules, &
119 local_particles=local_particles, &
120 molecules=molecules, &
121 molecule_kinds=molecule_kinds, &
122 particles=particles, &
123 gci=gci)
124 nparticle_kind = atomic_kinds%n_els
125 IF (PRESENT(compold)) THEN
126 IF (compold) THEN
127 CALL getold(gci, local_molecules, molecules%els, molecule_kinds%els, &
128 particles%els, cell)
129 END IF
130 END IF
131 has_pos = .false.
132 IF (.NOT. ASSOCIATED(my_pos)) THEN
133 has_pos = .true.
134 ALLOCATE (my_pos(3, particles%n_els))
135 my_pos = 0.0_dp
136 DO iparticle_kind = 1, nparticle_kind
137 nparticle_local = local_particles%n_el(iparticle_kind)
138 DO iparticle_local = 1, nparticle_local
139 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
140 my_pos(:, iparticle) = particles%els(iparticle)%r(:)
141 END DO
142 END DO
143 END IF
144 has_vel = .false.
145 IF (.NOT. ASSOCIATED(my_vel)) THEN
146 has_vel = .true.
147 ALLOCATE (my_vel(3, particles%n_els))
148 my_vel = 0.0_dp
149 DO iparticle_kind = 1, nparticle_kind
150 nparticle_local = local_particles%n_el(iparticle_kind)
151 DO iparticle_local = 1, nparticle_local
152 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
153 my_vel(:, iparticle) = particles%els(iparticle)%v(:)
154 END DO
155 END DO
156 END IF
157
158 CALL shake_control(gci=gci, local_molecules=local_molecules, &
159 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, &
160 particle_set=particles%els, pos=my_pos, vel=my_vel, dt=mydt, &
161 shake_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, &
162 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env, &
163 local_particles=local_particles)
164
165 ! Possibly reset the lagrange multipliers
166 IF (PRESENT(reset)) THEN
167 IF (reset) THEN
168 ! Reset Intramolecular constraints
169 DO i = 1, SIZE(molecules%els)
170 IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN
171 DO j = 1, SIZE(molecules%els(i)%lci%lcolv)
172 ! Reset langrange multiplier
173 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
174 END DO
175 END IF
176 IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN
177 DO j = 1, SIZE(molecules%els(i)%lci%lg3x3)
178 ! Reset langrange multiplier
179 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
180 END DO
181 END IF
182 IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN
183 DO j = 1, SIZE(molecules%els(i)%lci%lg4x6)
184 ! Reset langrange multiplier
185 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
186 END DO
187 END IF
188 END DO
189 ! Reset Intermolecular constraints
190 IF (ASSOCIATED(gci)) THEN
191 IF (ASSOCIATED(gci%lcolv)) THEN
192 DO j = 1, SIZE(gci%lcolv)
193 ! Reset langrange multiplier
194 gci%lcolv(j)%lambda = 0.0_dp
195 END DO
196 END IF
197 IF (ASSOCIATED(gci%lg3x3)) THEN
198 DO j = 1, SIZE(gci%lg3x3)
199 ! Reset langrange multiplier
200 gci%lg3x3(j)%lambda = 0.0_dp
201 END DO
202 END IF
203 IF (ASSOCIATED(gci%lg4x6)) THEN
204 DO j = 1, SIZE(gci%lg4x6)
205 ! Reset langrange multiplier
206 gci%lg4x6(j)%lambda = 0.0_dp
207 END DO
208 END IF
209 END IF
210 END IF
211 END IF
212
213 IF (has_pos) THEN
214 CALL update_particle_set(particles%els, force_env%para_env, pos=my_pos)
215 DEALLOCATE (my_pos)
216 END IF
217 IF (has_vel) THEN
218 CALL update_particle_set(particles%els, force_env%para_env, vel=my_vel)
219 DEALLOCATE (my_vel)
220 END IF
221 CALL timestop(handle)
222 END SUBROUTINE force_env_shake
223
224! **************************************************************************************************
225!> \brief perform rattle (enforcing of constraints on velocities)
226!> This routine can be easily adapted to performe rattle on whatever
227!> other vector different from forces..
228!> \param force_env the force env to shake
229!> \param dt the dt for shake (if you are not interested in the velocities
230!> it can be any positive number)
231!> \param shake_tol the tolerance for shake
232!> \param log_unit if >0 then some information on the shake is printed,
233!> defaults to -1
234!> \param lagrange_mult ...
235!> \param dump_lm ...
236!> \param vel ...
237!> \param reset ...
238!> \author tlaino
239! **************************************************************************************************
240 SUBROUTINE force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
241 vel, reset)
242
243 TYPE(force_env_type), POINTER :: force_env
244 REAL(kind=dp), INTENT(in), OPTIONAL :: dt
245 REAL(kind=dp), INTENT(in) :: shake_tol
246 INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult
247 LOGICAL, INTENT(IN), OPTIONAL :: dump_lm
248 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
249 OPTIONAL, TARGET :: vel
250 LOGICAL, INTENT(IN), OPTIONAL :: reset
251
252 CHARACTER(len=*), PARAMETER :: routinen = 'force_env_rattle'
253
254 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, &
255 my_log_unit, nparticle_kind, nparticle_local
256 LOGICAL :: has_vel, my_dump_lm
257 REAL(kind=dp) :: mydt
258 REAL(kind=dp), DIMENSION(:, :), POINTER :: my_vel
259 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
260 TYPE(cell_type), POINTER :: cell
261 TYPE(cp_subsys_type), POINTER :: subsys
262 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
263 TYPE(global_constraint_type), POINTER :: gci
264 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
265 TYPE(molecule_list_type), POINTER :: molecules
266 TYPE(particle_list_type), POINTER :: particles
267
268 CALL timeset(routinen, handle)
269 cpassert(ASSOCIATED(force_env))
270 cpassert(force_env%ref_count > 0)
271 my_log_unit = -1
272 IF (PRESENT(log_unit)) my_log_unit = log_unit
273 my_lagrange_mult = -1
274 IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
275 my_dump_lm = .false.
276 IF (PRESENT(dump_lm)) my_dump_lm = dump_lm
277 NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
278 my_vel)
279 IF (PRESENT(vel)) my_vel => vel
280 mydt = 0.1_dp
281 IF (PRESENT(dt)) mydt = dt
282 CALL force_env_get(force_env, subsys=subsys, cell=cell)
283 CALL cp_subsys_get(subsys, &
284 atomic_kinds=atomic_kinds, &
285 local_molecules=local_molecules, &
286 local_particles=local_particles, &
287 molecules=molecules, &
288 molecule_kinds=molecule_kinds, &
289 particles=particles, &
290 gci=gci)
291 nparticle_kind = atomic_kinds%n_els
292 has_vel = .false.
293 IF (.NOT. ASSOCIATED(my_vel)) THEN
294 has_vel = .true.
295 ALLOCATE (my_vel(3, particles%n_els))
296 my_vel = 0.0_dp
297 DO iparticle_kind = 1, nparticle_kind
298 nparticle_local = local_particles%n_el(iparticle_kind)
299 DO iparticle_local = 1, nparticle_local
300 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
301 my_vel(:, iparticle) = particles%els(iparticle)%v(:)
302 END DO
303 END DO
304 END IF
305
306 CALL rattle_control(gci=gci, local_molecules=local_molecules, &
307 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, &
308 particle_set=particles%els, vel=my_vel, dt=mydt, &
309 rattle_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, &
310 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env, &
311 local_particles=local_particles)
312
313 ! Possibly reset the lagrange multipliers
314 IF (PRESENT(reset)) THEN
315 IF (reset) THEN
316 ! Reset Intramolecular constraints
317 DO i = 1, SIZE(molecules%els)
318 IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN
319 DO j = 1, SIZE(molecules%els(i)%lci%lcolv)
320 ! Reset langrange multiplier
321 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
322 END DO
323 END IF
324 IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN
325 DO j = 1, SIZE(molecules%els(i)%lci%lg3x3)
326 ! Reset langrange multiplier
327 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
328 END DO
329 END IF
330 IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN
331 DO j = 1, SIZE(molecules%els(i)%lci%lg4x6)
332 ! Reset langrange multiplier
333 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
334 END DO
335 END IF
336 END DO
337 ! Reset Intermolecular constraints
338 IF (ASSOCIATED(gci)) THEN
339 IF (ASSOCIATED(gci%lcolv)) THEN
340 DO j = 1, SIZE(gci%lcolv)
341 ! Reset langrange multiplier
342 gci%lcolv(j)%lambda = 0.0_dp
343 END DO
344 END IF
345 IF (ASSOCIATED(gci%lg3x3)) THEN
346 DO j = 1, SIZE(gci%lg3x3)
347 ! Reset langrange multiplier
348 gci%lg3x3(j)%lambda = 0.0_dp
349 END DO
350 END IF
351 IF (ASSOCIATED(gci%lg4x6)) THEN
352 DO j = 1, SIZE(gci%lg4x6)
353 ! Reset langrange multiplier
354 gci%lg4x6(j)%lambda = 0.0_dp
355 END DO
356 END IF
357 END IF
358 END IF
359 END IF
360
361 IF (has_vel) THEN
362 CALL update_particle_set(particles%els, force_env%para_env, vel=my_vel)
363 END IF
364 DEALLOCATE (my_vel)
365 CALL timestop(handle)
366 END SUBROUTINE force_env_rattle
367
368! **************************************************************************************************
369!> \brief Rescale forces if requested
370!> \param force_env the force env to shake
371!> \author tlaino
372! **************************************************************************************************
373 SUBROUTINE rescale_forces(force_env)
374 TYPE(force_env_type), POINTER :: force_env
375
376 CHARACTER(len=*), PARAMETER :: routinen = 'rescale_forces'
377
378 INTEGER :: handle, iparticle
379 LOGICAL :: explicit
380 REAL(kind=dp) :: force(3), max_value, mod_force
381 TYPE(cp_subsys_type), POINTER :: subsys
382 TYPE(particle_list_type), POINTER :: particles
383 TYPE(section_vals_type), POINTER :: rescale_force_section
384
385 CALL timeset(routinen, handle)
386 cpassert(ASSOCIATED(force_env))
387 cpassert(force_env%ref_count > 0)
388 rescale_force_section => section_vals_get_subs_vals(force_env%force_env_section, "RESCALE_FORCES")
389 CALL section_vals_get(rescale_force_section, explicit=explicit)
390 IF (explicit) THEN
391 CALL section_vals_val_get(rescale_force_section, "MAX_FORCE", r_val=max_value)
392 CALL force_env_get(force_env, subsys=subsys)
393 CALL cp_subsys_get(subsys, particles=particles)
394 DO iparticle = 1, SIZE(particles%els)
395 force = particles%els(iparticle)%f(:)
396 mod_force = sqrt(dot_product(force, force))
397 IF ((mod_force > max_value) .AND. (mod_force /= 0.0_dp)) THEN
398 force = force/mod_force*max_value
399 particles%els(iparticle)%f(:) = force
400 END IF
401 END DO
402 END IF
403 CALL timestop(handle)
404 END SUBROUTINE rescale_forces
405
406! **************************************************************************************************
407!> \brief Write forces either to the screen or to a file.
408!> \param particles ...
409!> \param iw ...
410!> \param label ...
411!> \param ndigits ...
412!> \param unit_string ...
413!> \param total_force ...
414!> \param grand_total_force ...
415!> \param zero_force_core_shell_atom ...
416!> \author Ole Schuett
417! **************************************************************************************************
418 SUBROUTINE write_forces(particles, iw, label, ndigits, unit_string, total_force, &
419 grand_total_force, zero_force_core_shell_atom)
420
421 TYPE(particle_list_type), POINTER :: particles
422 INTEGER, INTENT(IN) :: iw
423 CHARACTER(LEN=*), INTENT(IN) :: label
424 INTEGER, INTENT(IN) :: ndigits
425 CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
426 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: total_force
427 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
428 OPTIONAL :: grand_total_force
429 LOGICAL, INTENT(IN), OPTIONAL :: zero_force_core_shell_atom
430
431 IF (iw == cp_logger_get_default_io_unit()) THEN
432 CALL write_forces_to_screen(particles, iw, label, ndigits, unit_string, total_force, &
433 grand_total_force, zero_force_core_shell_atom)
434 ELSE
435 CALL write_forces_to_file(particles, iw, label, ndigits, total_force, &
436 grand_total_force, zero_force_core_shell_atom)
437 END IF
438 END SUBROUTINE write_forces
439
440! **************************************************************************************************
441!> \brief Write forces to the screen.
442!>
443!> \param particles ...
444!> \param iw ...
445!> \param label ...
446!> \param ndigits ...
447!> \param unit_string ...
448!> \param total_force ...
449!> \param grand_total_force ...
450!> \param zero_force_core_shell_atom ...
451!> \author MK (06.09.2010)
452! **************************************************************************************************
453 SUBROUTINE write_forces_to_screen(particles, iw, label, ndigits, unit_string, total_force, &
454 grand_total_force, zero_force_core_shell_atom)
455
456 TYPE(particle_list_type), POINTER :: particles
457 INTEGER, INTENT(IN) :: iw
458 CHARACTER(LEN=*), INTENT(IN) :: label
459 INTEGER, INTENT(IN) :: ndigits
460 CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
461 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: total_force
462 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
463 OPTIONAL :: grand_total_force
464 LOGICAL, INTENT(IN), OPTIONAL :: zero_force_core_shell_atom
465
466 CHARACTER(LEN=18) :: fmtstr4
467 CHARACTER(LEN=29) :: fmtstr3
468 CHARACTER(LEN=38) :: fmtstr2
469 CHARACTER(LEN=49) :: fmtstr1
470 CHARACTER(LEN=7) :: tag
471 CHARACTER(LEN=LEN_TRIM(label)) :: lc_label
472 INTEGER :: i, iparticle, n
473 LOGICAL :: zero_force
474 REAL(kind=dp) :: fconv
475 REAL(kind=dp), DIMENSION(3) :: f
476
477 IF (iw > 0) THEN
478 cpassert(ASSOCIATED(particles))
479 tag = "FORCES|"
480 lc_label = trim(label)
481 CALL lowercase(lc_label)
482 n = min(max(1, ndigits), 20)
483 fmtstr1 = "(/,T2,A,1X,A,/,T2,A,3X,A,T20,A3,2( X,A3), X,A3)"
484 WRITE (unit=fmtstr1(35:36), fmt="(I2)") n + 5
485 WRITE (unit=fmtstr1(43:44), fmt="(I2)") n + 6
486 fmtstr2 = "(T2,A,I7,T16,3(1X,ES . ),2X,ES . )"
487 WRITE (unit=fmtstr2(21:22), fmt="(I2)") n + 7
488 WRITE (unit=fmtstr2(24:25), fmt="(I2)") n
489 WRITE (unit=fmtstr2(33:34), fmt="(I2)") n + 7
490 WRITE (unit=fmtstr2(36:37), fmt="(I2)") n
491 fmtstr3 = "(T2,A,T16,3(1X,ES . ))"
492 WRITE (unit=fmtstr3(18:19), fmt="(I2)") n + 7
493 WRITE (unit=fmtstr3(21:22), fmt="(I2)") n
494 fmtstr4 = "(T2,A,T ,ES . )"
495 WRITE (unit=fmtstr4(8:9), fmt="(I2)") 3*(n + 8) + 18
496 WRITE (unit=fmtstr4(13:14), fmt="(I2)") n + 7
497 WRITE (unit=fmtstr4(16:17), fmt="(I2)") n
498 IF (PRESENT(zero_force_core_shell_atom)) THEN
499 zero_force = zero_force_core_shell_atom
500 ELSE
501 zero_force = .false.
502 END IF
503 fconv = cp_unit_from_cp2k(1.0_dp, trim(unit_string))
504 WRITE (unit=iw, fmt=fmtstr1) &
505 tag, label//" forces ["//trim(adjustl(unit_string))//"]", &
506 tag, "Atom", " x ", " y ", " z ", "|f|"
507 total_force(1:3) = 0.0_dp
508 DO iparticle = 1, particles%n_els
509 IF (particles%els(iparticle)%atom_index /= 0) THEN
510 i = particles%els(iparticle)%atom_index
511 ELSE
512 i = iparticle
513 END IF
514 IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0)) THEN
515 f(1:3) = 0.0_dp
516 ELSE
517 f(1:3) = particles%els(iparticle)%f(1:3)*fconv
518 END IF
519 WRITE (unit=iw, fmt=fmtstr2) &
520 tag, i, f(1:3), sqrt(sum(f(1:3)**2))
521 total_force(1:3) = total_force(1:3) + f(1:3)
522 END DO
523 WRITE (unit=iw, fmt=fmtstr3) &
524 tag//" Sum", total_force(1:3)
525 WRITE (unit=iw, fmt=fmtstr4) &
526 tag//" Total "//trim(adjustl(lc_label))//" force", &
527 sqrt(sum(total_force(1:3)**2))
528 END IF
529
530 IF (PRESENT(grand_total_force)) THEN
531 grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
532 WRITE (unit=iw, fmt="(A)") ""
533 WRITE (unit=iw, fmt=fmtstr4) &
534 tag//" Grand total force ["//trim(adjustl(unit_string))//"]", &
535 sqrt(sum(grand_total_force(1:3)**2))
536 END IF
537
538 END SUBROUTINE write_forces_to_screen
539
540! **************************************************************************************************
541!> \brief Write forces to a file using our stable legacy format. Please don't change the format.
542!>
543!> \param particles ...
544!> \param iw ...
545!> \param label ...
546!> \param ndigits ...
547!> \param total_force ...
548!> \param grand_total_force ...
549!> \param zero_force_core_shell_atom ...
550!> \author MK (06.09.2010)
551! **************************************************************************************************
552 SUBROUTINE write_forces_to_file(particles, iw, label, ndigits, total_force, &
553 grand_total_force, zero_force_core_shell_atom)
554
555 TYPE(particle_list_type), POINTER :: particles
556 INTEGER, INTENT(IN) :: iw
557 CHARACTER(LEN=*), INTENT(IN) :: label
558 INTEGER, INTENT(IN) :: ndigits
559 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: total_force
560 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
561 OPTIONAL :: grand_total_force
562 LOGICAL, INTENT(IN), OPTIONAL :: zero_force_core_shell_atom
563
564 CHARACTER(LEN=23) :: fmtstr3
565 CHARACTER(LEN=36) :: fmtstr2
566 CHARACTER(LEN=46) :: fmtstr1
567 CHARACTER(LEN=LEN_TRIM(label)) :: uc_label
568 INTEGER :: i, ikind, iparticle, n
569 LOGICAL :: zero_force
570 REAL(kind=dp), DIMENSION(3) :: f
571
572 IF (iw > 0) THEN
573 cpassert(ASSOCIATED(particles))
574 n = min(max(1, ndigits), 20)
575 fmtstr1 = "(/,T2,A,/,/,T2,A,T11,A,T18,A,T35,A1,2( X,A1))"
576 WRITE (unit=fmtstr1(39:40), fmt="(I2)") n + 6
577 fmtstr2 = "(T2,I6,1X,I6,T21,A,T28,3(1X,F . ))"
578 WRITE (unit=fmtstr2(33:34), fmt="(I2)") n
579 WRITE (unit=fmtstr2(30:31), fmt="(I2)") n + 6
580 fmtstr3 = "(T2,A,T28,4(1X,F . ))"
581 WRITE (unit=fmtstr3(20:21), fmt="(I2)") n
582 WRITE (unit=fmtstr3(17:18), fmt="(I2)") n + 6
583 IF (PRESENT(zero_force_core_shell_atom)) THEN
584 zero_force = zero_force_core_shell_atom
585 ELSE
586 zero_force = .false.
587 END IF
588 uc_label = trim(label)
589 CALL uppercase(uc_label)
590 WRITE (unit=iw, fmt=fmtstr1) &
591 uc_label//" FORCES in [a.u.]", "# Atom", "Kind", "Element", "X", "Y", "Z"
592 total_force(1:3) = 0.0_dp
593 DO iparticle = 1, particles%n_els
594 ikind = particles%els(iparticle)%atomic_kind%kind_number
595 IF (particles%els(iparticle)%atom_index /= 0) THEN
596 i = particles%els(iparticle)%atom_index
597 ELSE
598 i = iparticle
599 END IF
600 IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0)) THEN
601 f(1:3) = 0.0_dp
602 ELSE
603 f(1:3) = particles%els(iparticle)%f(1:3)
604 END IF
605 WRITE (unit=iw, fmt=fmtstr2) &
606 i, ikind, particles%els(iparticle)%atomic_kind%element_symbol, f(1:3)
607 total_force(1:3) = total_force(1:3) + f(1:3)
608 END DO
609 WRITE (unit=iw, fmt=fmtstr3) &
610 "SUM OF "//uc_label//" FORCES", total_force(1:3), sqrt(sum(total_force(:)**2))
611 END IF
612
613 IF (PRESENT(grand_total_force)) THEN
614 grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
615 WRITE (unit=iw, fmt="(A)") ""
616 WRITE (unit=iw, fmt=fmtstr3) &
617 "GRAND TOTAL FORCE", grand_total_force(1:3), sqrt(sum(grand_total_force(:)**2))
618 END IF
619
620 END SUBROUTINE write_forces_to_file
621
622! **************************************************************************************************
623!> \brief Write the atomic coordinates, types, and energies
624!> \param iounit ...
625!> \param particles ...
626!> \param atener ...
627!> \param label ...
628!> \date 05.06.2023
629!> \author JGH
630!> \version 1.0
631! **************************************************************************************************
632 SUBROUTINE write_atener(iounit, particles, atener, label)
633
634 INTEGER, INTENT(IN) :: iounit
635 TYPE(particle_list_type), POINTER :: particles
636 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: atener
637 CHARACTER(LEN=*), INTENT(IN) :: label
638
639 INTEGER :: i, natom
640
641 IF (iounit > 0) THEN
642 WRITE (unit=iounit, fmt="(/,T2,A)") trim(label)
643 WRITE (unit=iounit, fmt="(T4,A,T30,A,T42,A,T54,A,T69,A)") &
644 "Atom Element", "X", "Y", "Z", "Energy[a.u.]"
645 natom = particles%n_els
646 DO i = 1, natom
647 WRITE (iounit, "(I6,T12,A2,T24,3F12.6,F21.10)") i, &
648 trim(adjustl(particles%els(i)%atomic_kind%element_symbol)), &
649 particles%els(i)%r(1:3)*angstrom, atener(i)
650 END DO
651 WRITE (iounit, "(A)") ""
652 END IF
653
654 END SUBROUTINE write_atener
655
656END MODULE force_env_utils
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
constraint_util
Contains routines useful for the application of constraints during MD.
Definition constraint_util.F:13
constraint_util::getold
subroutine, public getold(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, cell)
saves all of the old variables
Definition constraint_util.F:63
constraint::rattle_control
subroutine, public rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, vel, dt, rattle_tol, log_unit, lagrange_mult, dump_lm, cell, group, local_particles)
...
Definition constraint.F:229
constraint::shake_control
subroutine, public shake_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, pos, vel, dt, shake_tol, log_unit, lagrange_mult, dump_lm, cell, group, local_particles)
...
Definition constraint.F:101
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1178
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
force_env_utils
Util force_env module.
Definition force_env_utils.F:12
force_env_utils::write_forces
subroutine, public write_forces(particles, iw, label, ndigits, unit_string, total_force, grand_total_force, zero_force_core_shell_atom)
Write forces either to the screen or to a file.
Definition force_env_utils.F:420
force_env_utils::force_env_shake
subroutine, public force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, pos, vel, compold, reset)
perform shake (enforcing of constraints)
Definition force_env_utils.F:74
force_env_utils::write_atener
subroutine, public write_atener(iounit, particles, atener, label)
Write the atomic coordinates, types, and energies.
Definition force_env_utils.F:633
force_env_utils::rescale_forces
subroutine, public rescale_forces(force_env)
Rescale forces if requested.
Definition force_env_utils.F:374
force_env_utils::force_env_rattle
subroutine, public force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, vel, reset)
perform rattle (enforcing of constraints on velocities) This routine can be easily adapted to perform...
Definition force_env_utils.F:242
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_list_types
represent a simple array based list of the given type
Definition molecule_list_types.F:15
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
particle_types::update_particle_set
subroutine, public update_particle_set(particle_set, int_group, pos, vel, for, add)
...
Definition particle_types.F:106
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::lowercase
elemental subroutine, public lowercase(string)
Convert all upper case characters in a string to lower case.
Definition string_utilities.F:3341
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::global_constraint_type
Definition molecule_types.F:93
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46