Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Util force_env module
10 : !> \author Teodoro Laino [tlaino] - 02.2011
11 : ! **************************************************************************************************
12 : MODULE 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 :
54 : CONTAINS
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 48 : SUBROUTINE force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
73 24 : 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 24 : 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 24 : CALL timeset(routineN, handle)
101 24 : CPASSERT(ASSOCIATED(force_env))
102 24 : CPASSERT(force_env%ref_count > 0)
103 24 : my_log_unit = -1
104 24 : IF (PRESENT(log_unit)) my_log_unit = log_unit
105 24 : my_lagrange_mult = -1
106 24 : IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
107 24 : my_dump_lm = .FALSE.
108 24 : IF (PRESENT(dump_lm)) my_dump_lm = dump_lm
109 24 : NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
110 24 : my_pos, my_vel, gci)
111 24 : IF (PRESENT(pos)) my_pos => pos
112 24 : IF (PRESENT(vel)) my_vel => vel
113 24 : mydt = 0.1_dp
114 24 : IF (PRESENT(dt)) mydt = dt
115 24 : 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 24 : gci=gci)
124 24 : nparticle_kind = atomic_kinds%n_els
125 24 : IF (PRESENT(compold)) THEN
126 24 : IF (compold) THEN
127 : CALL getold(gci, local_molecules, molecules%els, molecule_kinds%els, &
128 24 : particles%els, cell)
129 : END IF
130 : END IF
131 24 : has_pos = .FALSE.
132 24 : IF (.NOT. ASSOCIATED(my_pos)) THEN
133 24 : has_pos = .TRUE.
134 72 : ALLOCATE (my_pos(3, particles%n_els))
135 7504 : my_pos = 0.0_dp
136 112 : DO iparticle_kind = 1, nparticle_kind
137 88 : nparticle_local = local_particles%n_el(iparticle_kind)
138 1047 : DO iparticle_local = 1, nparticle_local
139 935 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
140 3828 : my_pos(:, iparticle) = particles%els(iparticle)%r(:)
141 : END DO
142 : END DO
143 : END IF
144 24 : has_vel = .FALSE.
145 24 : IF (.NOT. ASSOCIATED(my_vel)) THEN
146 24 : has_vel = .TRUE.
147 72 : ALLOCATE (my_vel(3, particles%n_els))
148 7504 : my_vel = 0.0_dp
149 112 : DO iparticle_kind = 1, nparticle_kind
150 88 : nparticle_local = local_particles%n_el(iparticle_kind)
151 1047 : DO iparticle_local = 1, nparticle_local
152 935 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
153 3828 : 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 24 : local_particles=local_particles)
164 :
165 : ! Possibly reset the lagrange multipliers
166 24 : IF (PRESENT(reset)) THEN
167 0 : IF (reset) THEN
168 : ! Reset Intramolecular constraints
169 0 : DO i = 1, SIZE(molecules%els)
170 0 : IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN
171 0 : DO j = 1, SIZE(molecules%els(i)%lci%lcolv)
172 : ! Reset langrange multiplier
173 0 : molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
174 : END DO
175 : END IF
176 0 : IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN
177 0 : DO j = 1, SIZE(molecules%els(i)%lci%lg3x3)
178 : ! Reset langrange multiplier
179 0 : molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
180 : END DO
181 : END IF
182 0 : IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN
183 0 : DO j = 1, SIZE(molecules%els(i)%lci%lg4x6)
184 : ! Reset langrange multiplier
185 0 : molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
186 : END DO
187 : END IF
188 : END DO
189 : ! Reset Intermolecular constraints
190 0 : IF (ASSOCIATED(gci)) THEN
191 0 : IF (ASSOCIATED(gci%lcolv)) THEN
192 0 : DO j = 1, SIZE(gci%lcolv)
193 : ! Reset langrange multiplier
194 0 : gci%lcolv(j)%lambda = 0.0_dp
195 : END DO
196 : END IF
197 0 : IF (ASSOCIATED(gci%lg3x3)) THEN
198 0 : DO j = 1, SIZE(gci%lg3x3)
199 : ! Reset langrange multiplier
200 0 : gci%lg3x3(j)%lambda = 0.0_dp
201 : END DO
202 : END IF
203 0 : IF (ASSOCIATED(gci%lg4x6)) THEN
204 0 : DO j = 1, SIZE(gci%lg4x6)
205 : ! Reset langrange multiplier
206 0 : gci%lg4x6(j)%lambda = 0.0_dp
207 : END DO
208 : END IF
209 : END IF
210 : END IF
211 : END IF
212 :
213 24 : IF (has_pos) THEN
214 24 : CALL update_particle_set(particles%els, force_env%para_env, pos=my_pos)
215 24 : DEALLOCATE (my_pos)
216 : END IF
217 24 : IF (has_vel) THEN
218 24 : CALL update_particle_set(particles%els, force_env%para_env, vel=my_vel)
219 24 : DEALLOCATE (my_vel)
220 : END IF
221 24 : CALL timestop(handle)
222 24 : 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 48 : SUBROUTINE force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
241 24 : 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 24 : 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 24 : CALL timeset(routineN, handle)
269 24 : CPASSERT(ASSOCIATED(force_env))
270 24 : CPASSERT(force_env%ref_count > 0)
271 24 : my_log_unit = -1
272 24 : IF (PRESENT(log_unit)) my_log_unit = log_unit
273 24 : my_lagrange_mult = -1
274 24 : IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
275 24 : my_dump_lm = .FALSE.
276 24 : IF (PRESENT(dump_lm)) my_dump_lm = dump_lm
277 24 : NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
278 24 : my_vel)
279 24 : IF (PRESENT(vel)) my_vel => vel
280 24 : mydt = 0.1_dp
281 24 : IF (PRESENT(dt)) mydt = dt
282 24 : 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 24 : gci=gci)
291 24 : nparticle_kind = atomic_kinds%n_els
292 24 : has_vel = .FALSE.
293 24 : IF (.NOT. ASSOCIATED(my_vel)) THEN
294 24 : has_vel = .TRUE.
295 72 : ALLOCATE (my_vel(3, particles%n_els))
296 7504 : my_vel = 0.0_dp
297 112 : DO iparticle_kind = 1, nparticle_kind
298 88 : nparticle_local = local_particles%n_el(iparticle_kind)
299 1047 : DO iparticle_local = 1, nparticle_local
300 935 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
301 3828 : 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 24 : local_particles=local_particles)
312 :
313 : ! Possibly reset the lagrange multipliers
314 24 : IF (PRESENT(reset)) THEN
315 24 : IF (reset) THEN
316 : ! Reset Intramolecular constraints
317 500 : DO i = 1, SIZE(molecules%els)
318 476 : IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN
319 28 : DO j = 1, SIZE(molecules%els(i)%lci%lcolv)
320 : ! Reset langrange multiplier
321 28 : molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
322 : END DO
323 : END IF
324 476 : IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN
325 890 : DO j = 1, SIZE(molecules%els(i)%lci%lg3x3)
326 : ! Reset langrange multiplier
327 2216 : molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
328 : END DO
329 : END IF
330 500 : IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN
331 12 : DO j = 1, SIZE(molecules%els(i)%lci%lg4x6)
332 : ! Reset langrange multiplier
333 36 : molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
334 : END DO
335 : END IF
336 : END DO
337 : ! Reset Intermolecular constraints
338 24 : IF (ASSOCIATED(gci)) THEN
339 24 : IF (ASSOCIATED(gci%lcolv)) THEN
340 0 : DO j = 1, SIZE(gci%lcolv)
341 : ! Reset langrange multiplier
342 0 : gci%lcolv(j)%lambda = 0.0_dp
343 : END DO
344 : END IF
345 24 : IF (ASSOCIATED(gci%lg3x3)) THEN
346 0 : DO j = 1, SIZE(gci%lg3x3)
347 : ! Reset langrange multiplier
348 0 : gci%lg3x3(j)%lambda = 0.0_dp
349 : END DO
350 : END IF
351 24 : IF (ASSOCIATED(gci%lg4x6)) THEN
352 0 : DO j = 1, SIZE(gci%lg4x6)
353 : ! Reset langrange multiplier
354 0 : gci%lg4x6(j)%lambda = 0.0_dp
355 : END DO
356 : END IF
357 : END IF
358 : END IF
359 : END IF
360 :
361 24 : IF (has_vel) THEN
362 24 : CALL update_particle_set(particles%els, force_env%para_env, vel=my_vel)
363 : END IF
364 24 : DEALLOCATE (my_vel)
365 24 : CALL timestop(handle)
366 24 : 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 195142 : 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 97571 : CALL timeset(routineN, handle)
386 97571 : CPASSERT(ASSOCIATED(force_env))
387 97571 : CPASSERT(force_env%ref_count > 0)
388 97571 : rescale_force_section => section_vals_get_subs_vals(force_env%force_env_section, "RESCALE_FORCES")
389 97571 : CALL section_vals_get(rescale_force_section, explicit=explicit)
390 97571 : IF (explicit) THEN
391 224 : CALL section_vals_val_get(rescale_force_section, "MAX_FORCE", r_val=max_value)
392 224 : CALL force_env_get(force_env, subsys=subsys)
393 224 : CALL cp_subsys_get(subsys, particles=particles)
394 36262 : DO iparticle = 1, SIZE(particles%els)
395 144152 : force = particles%els(iparticle)%f(:)
396 144152 : mod_force = SQRT(DOT_PRODUCT(force, force))
397 36262 : IF ((mod_force > max_value) .AND. (mod_force /= 0.0_dp)) THEN
398 10208 : force = force/mod_force*max_value
399 10208 : particles%els(iparticle)%f(:) = force
400 : END IF
401 : END DO
402 : END IF
403 97571 : CALL timestop(handle)
404 97571 : 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 1835 : 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 1835 : IF (iw == cp_logger_get_default_io_unit()) THEN
432 : CALL write_forces_to_screen(particles, iw, label, ndigits, unit_string, total_force, &
433 1830 : grand_total_force, zero_force_core_shell_atom)
434 : ELSE
435 : CALL write_forces_to_file(particles, iw, label, ndigits, total_force, &
436 5 : grand_total_force, zero_force_core_shell_atom)
437 : END IF
438 1835 : 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 1830 : 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 1830 : 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 1830 : IF (iw > 0) THEN
478 1830 : CPASSERT(ASSOCIATED(particles))
479 1830 : tag = "FORCES|"
480 1830 : lc_label = TRIM(label)
481 1830 : CALL lowercase(lc_label)
482 1830 : n = MIN(MAX(1, ndigits), 20)
483 1830 : fmtstr1 = "(/,T2,A,1X,A,/,T2,A,3X,A,T20,A3,2( X,A3), X,A3)"
484 1830 : WRITE (UNIT=fmtstr1(35:36), FMT="(I2)") n + 5
485 1830 : WRITE (UNIT=fmtstr1(43:44), FMT="(I2)") n + 6
486 1830 : fmtstr2 = "(T2,A,I7,T16,3(1X,ES . ),2X,ES . )"
487 1830 : WRITE (UNIT=fmtstr2(21:22), FMT="(I2)") n + 7
488 1830 : WRITE (UNIT=fmtstr2(24:25), FMT="(I2)") n
489 1830 : WRITE (UNIT=fmtstr2(33:34), FMT="(I2)") n + 7
490 1830 : WRITE (UNIT=fmtstr2(36:37), FMT="(I2)") n
491 1830 : fmtstr3 = "(T2,A,T16,3(1X,ES . ))"
492 1830 : WRITE (UNIT=fmtstr3(18:19), FMT="(I2)") n + 7
493 1830 : WRITE (UNIT=fmtstr3(21:22), FMT="(I2)") n
494 1830 : fmtstr4 = "(T2,A,T ,ES . )"
495 1830 : WRITE (UNIT=fmtstr4(8:9), FMT="(I2)") 3*(n + 8) + 18
496 1830 : WRITE (UNIT=fmtstr4(13:14), FMT="(I2)") n + 7
497 1830 : WRITE (UNIT=fmtstr4(16:17), FMT="(I2)") n
498 1830 : IF (PRESENT(zero_force_core_shell_atom)) THEN
499 495 : zero_force = zero_force_core_shell_atom
500 : ELSE
501 : zero_force = .FALSE.
502 : END IF
503 1830 : fconv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_string))
504 : WRITE (UNIT=iw, FMT=fmtstr1) &
505 1830 : tag, label//" forces ["//TRIM(ADJUSTL(unit_string))//"]", &
506 3660 : tag, "Atom", " x ", " y ", " z ", "|f|"
507 1830 : total_force(1:3) = 0.0_dp
508 188619 : DO iparticle = 1, particles%n_els
509 186789 : IF (particles%els(iparticle)%atom_index /= 0) THEN
510 186789 : i = particles%els(iparticle)%atom_index
511 : ELSE
512 0 : i = iparticle
513 : END IF
514 186789 : IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0)) THEN
515 49172 : f(1:3) = 0.0_dp
516 : ELSE
517 550468 : f(1:3) = particles%els(iparticle)%f(1:3)*fconv
518 : END IF
519 : WRITE (UNIT=iw, FMT=fmtstr2) &
520 933945 : tag, i, f(1:3), SQRT(SUM(f(1:3)**2))
521 748986 : total_force(1:3) = total_force(1:3) + f(1:3)
522 : END DO
523 : WRITE (UNIT=iw, FMT=fmtstr3) &
524 1830 : tag//" Sum", total_force(1:3)
525 : WRITE (UNIT=iw, FMT=fmtstr4) &
526 1830 : tag//" Total "//TRIM(ADJUSTL(lc_label))//" force", &
527 9150 : SQRT(SUM(total_force(1:3)**2))
528 : END IF
529 :
530 1830 : IF (PRESENT(grand_total_force)) THEN
531 660 : grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
532 165 : WRITE (UNIT=iw, FMT="(A)") ""
533 : WRITE (UNIT=iw, FMT=fmtstr4) &
534 165 : tag//" Grand total force ["//TRIM(ADJUSTL(unit_string))//"]", &
535 825 : SQRT(SUM(grand_total_force(1:3)**2))
536 : END IF
537 :
538 1830 : 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 5 : 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 5 : 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 5 : IF (iw > 0) THEN
573 5 : CPASSERT(ASSOCIATED(particles))
574 5 : n = MIN(MAX(1, ndigits), 20)
575 5 : fmtstr1 = "(/,T2,A,/,/,T2,A,T11,A,T18,A,T35,A1,2( X,A1))"
576 5 : WRITE (UNIT=fmtstr1(39:40), FMT="(I2)") n + 6
577 5 : fmtstr2 = "(T2,I6,1X,I6,T21,A,T28,3(1X,F . ))"
578 5 : WRITE (UNIT=fmtstr2(33:34), FMT="(I2)") n
579 5 : WRITE (UNIT=fmtstr2(30:31), FMT="(I2)") n + 6
580 5 : fmtstr3 = "(T2,A,T28,4(1X,F . ))"
581 5 : WRITE (UNIT=fmtstr3(20:21), FMT="(I2)") n
582 5 : WRITE (UNIT=fmtstr3(17:18), FMT="(I2)") n + 6
583 5 : IF (PRESENT(zero_force_core_shell_atom)) THEN
584 0 : zero_force = zero_force_core_shell_atom
585 : ELSE
586 : zero_force = .FALSE.
587 : END IF
588 5 : uc_label = TRIM(label)
589 5 : CALL uppercase(uc_label)
590 : WRITE (UNIT=iw, FMT=fmtstr1) &
591 5 : uc_label//" FORCES in [a.u.]", "# Atom", "Kind", "Element", "X", "Y", "Z"
592 5 : total_force(1:3) = 0.0_dp
593 110 : DO iparticle = 1, particles%n_els
594 105 : ikind = particles%els(iparticle)%atomic_kind%kind_number
595 105 : IF (particles%els(iparticle)%atom_index /= 0) THEN
596 105 : i = particles%els(iparticle)%atom_index
597 : ELSE
598 0 : i = iparticle
599 : END IF
600 105 : IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0)) THEN
601 0 : f(1:3) = 0.0_dp
602 : ELSE
603 420 : f(1:3) = particles%els(iparticle)%f(1:3)
604 : END IF
605 : WRITE (UNIT=iw, FMT=fmtstr2) &
606 105 : i, ikind, particles%els(iparticle)%atomic_kind%element_symbol, f(1:3)
607 425 : total_force(1:3) = total_force(1:3) + f(1:3)
608 : END DO
609 : WRITE (UNIT=iw, FMT=fmtstr3) &
610 25 : "SUM OF "//uc_label//" FORCES", total_force(1:3), SQRT(SUM(total_force(:)**2))
611 : END IF
612 :
613 5 : IF (PRESENT(grand_total_force)) THEN
614 0 : grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
615 0 : WRITE (UNIT=iw, FMT="(A)") ""
616 : WRITE (UNIT=iw, FMT=fmtstr3) &
617 0 : "GRAND TOTAL FORCE", grand_total_force(1:3), SQRT(SUM(grand_total_force(:)**2))
618 : END IF
619 :
620 5 : 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 489 : 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 489 : IF (iounit > 0) THEN
642 489 : WRITE (UNIT=iounit, FMT="(/,T2,A)") TRIM(label)
643 : WRITE (UNIT=iounit, FMT="(T4,A,T30,A,T42,A,T54,A,T69,A)") &
644 489 : "Atom Element", "X", "Y", "Z", "Energy[a.u.]"
645 489 : natom = particles%n_els
646 17788 : DO i = 1, natom
647 17299 : WRITE (iounit, "(I6,T12,A2,T24,3F12.6,F21.10)") i, &
648 17299 : TRIM(ADJUSTL(particles%els(i)%atomic_kind%element_symbol)), &
649 86984 : particles%els(i)%r(1:3)*angstrom, atener(i)
650 : END DO
651 489 : WRITE (iounit, "(A)") ""
652 : END IF
653 :
654 489 : END SUBROUTINE write_atener
655 :
656 : END MODULE force_env_utils
|
