Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_subsys_types, ONLY: cp_subsys_get,&
20 : cp_subsys_type
21 : USE cp_units, ONLY: cp_unit_from_cp2k
22 : USE distribution_1d_types, ONLY: distribution_1d_type
23 : USE force_env_types, ONLY: force_env_get,&
24 : force_env_type
25 : USE input_section_types, ONLY: section_vals_get,&
26 : section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: default_string_length,&
30 : dp
31 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
32 : USE molecule_list_types, ONLY: molecule_list_type
33 : USE molecule_types, ONLY: global_constraint_type
34 : USE particle_list_types, ONLY: particle_list_type
35 : USE particle_types, ONLY: update_particle_set
36 : USE physcon, ONLY: angstrom
37 : USE string_utilities, ONLY: lowercase
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_utils'
45 :
46 : PUBLIC :: force_env_shake, &
47 : force_env_rattle, &
48 : rescale_forces, &
49 : write_atener, &
50 : write_forces
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief perform shake (enforcing of constraints)
56 : !> \param force_env the force env to shake
57 : !> \param dt the dt for shake (if you are not interested in the velocities
58 : !> it can be any positive number)
59 : !> \param shake_tol the tolerance for shake
60 : !> \param log_unit if >0 then some information on the shake is printed,
61 : !> defaults to -1
62 : !> \param lagrange_mult ...
63 : !> \param dump_lm ...
64 : !> \param pos ...
65 : !> \param vel ...
66 : !> \param compold ...
67 : !> \param reset ...
68 : !> \author fawzi
69 : ! **************************************************************************************************
70 48 : SUBROUTINE force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
71 24 : pos, vel, compold, reset)
72 :
73 : TYPE(force_env_type), POINTER :: force_env
74 : REAL(kind=dp), INTENT(IN), OPTIONAL :: dt
75 : REAL(kind=dp), INTENT(IN) :: shake_tol
76 : INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult
77 : LOGICAL, INTENT(IN), OPTIONAL :: dump_lm
78 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
79 : OPTIONAL, TARGET :: pos, vel
80 : LOGICAL, INTENT(IN), OPTIONAL :: compold, reset
81 :
82 : CHARACTER(len=*), PARAMETER :: routineN = 'force_env_shake'
83 :
84 : INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, &
85 : my_log_unit, nparticle_kind, nparticle_local
86 : LOGICAL :: has_pos, has_vel, my_dump_lm
87 : REAL(KIND=dp) :: mydt
88 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_pos, my_vel
89 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
90 : TYPE(cell_type), POINTER :: cell
91 : TYPE(cp_subsys_type), POINTER :: subsys
92 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
93 : TYPE(global_constraint_type), POINTER :: gci
94 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
95 : TYPE(molecule_list_type), POINTER :: molecules
96 : TYPE(particle_list_type), POINTER :: particles
97 :
98 24 : CALL timeset(routineN, handle)
99 24 : CPASSERT(ASSOCIATED(force_env))
100 24 : CPASSERT(force_env%ref_count > 0)
101 24 : my_log_unit = -1
102 24 : IF (PRESENT(log_unit)) my_log_unit = log_unit
103 24 : my_lagrange_mult = -1
104 24 : IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
105 24 : my_dump_lm = .FALSE.
106 24 : IF (PRESENT(dump_lm)) my_dump_lm = dump_lm
107 24 : NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
108 24 : my_pos, my_vel, gci)
109 24 : IF (PRESENT(pos)) my_pos => pos
110 24 : IF (PRESENT(vel)) my_vel => vel
111 24 : mydt = 0.1_dp
112 24 : IF (PRESENT(dt)) mydt = dt
113 24 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
114 : CALL cp_subsys_get(subsys, &
115 : atomic_kinds=atomic_kinds, &
116 : local_molecules=local_molecules, &
117 : local_particles=local_particles, &
118 : molecules=molecules, &
119 : molecule_kinds=molecule_kinds, &
120 : particles=particles, &
121 24 : gci=gci)
122 24 : nparticle_kind = atomic_kinds%n_els
123 24 : IF (PRESENT(compold)) THEN
124 24 : IF (compold) THEN
125 : CALL getold(gci, local_molecules, molecules%els, molecule_kinds%els, &
126 24 : particles%els, cell)
127 : END IF
128 : END IF
129 24 : has_pos = .FALSE.
130 24 : IF (.NOT. ASSOCIATED(my_pos)) THEN
131 24 : has_pos = .TRUE.
132 72 : ALLOCATE (my_pos(3, particles%n_els))
133 7504 : my_pos = 0.0_dp
134 112 : DO iparticle_kind = 1, nparticle_kind
135 88 : nparticle_local = local_particles%n_el(iparticle_kind)
136 1047 : DO iparticle_local = 1, nparticle_local
137 935 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
138 3828 : my_pos(:, iparticle) = particles%els(iparticle)%r(:)
139 : END DO
140 : END DO
141 : END IF
142 24 : has_vel = .FALSE.
143 24 : IF (.NOT. ASSOCIATED(my_vel)) THEN
144 24 : has_vel = .TRUE.
145 72 : ALLOCATE (my_vel(3, particles%n_els))
146 7504 : my_vel = 0.0_dp
147 112 : DO iparticle_kind = 1, nparticle_kind
148 88 : nparticle_local = local_particles%n_el(iparticle_kind)
149 1047 : DO iparticle_local = 1, nparticle_local
150 935 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
151 3828 : my_vel(:, iparticle) = particles%els(iparticle)%v(:)
152 : END DO
153 : END DO
154 : END IF
155 :
156 : CALL shake_control(gci=gci, local_molecules=local_molecules, &
157 : molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, &
158 : particle_set=particles%els, pos=my_pos, vel=my_vel, dt=mydt, &
159 : shake_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, &
160 : dump_lm=my_dump_lm, cell=cell, group=force_env%para_env, &
161 24 : local_particles=local_particles)
162 :
163 : ! Possibly reset the lagrange multipliers
164 24 : IF (PRESENT(reset)) THEN
165 0 : IF (reset) THEN
166 : ! Reset Intramolecular constraints
167 0 : DO i = 1, SIZE(molecules%els)
168 0 : IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN
169 0 : DO j = 1, SIZE(molecules%els(i)%lci%lcolv)
170 : ! Reset langrange multiplier
171 0 : molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
172 : END DO
173 : END IF
174 0 : IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN
175 0 : DO j = 1, SIZE(molecules%els(i)%lci%lg3x3)
176 : ! Reset langrange multiplier
177 0 : molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
178 : END DO
179 : END IF
180 0 : IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN
181 0 : DO j = 1, SIZE(molecules%els(i)%lci%lg4x6)
182 : ! Reset langrange multiplier
183 0 : molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
184 : END DO
185 : END IF
186 : END DO
187 : ! Reset Intermolecular constraints
188 0 : IF (ASSOCIATED(gci)) THEN
189 0 : IF (ASSOCIATED(gci%lcolv)) THEN
190 0 : DO j = 1, SIZE(gci%lcolv)
191 : ! Reset langrange multiplier
192 0 : gci%lcolv(j)%lambda = 0.0_dp
193 : END DO
194 : END IF
195 0 : IF (ASSOCIATED(gci%lg3x3)) THEN
196 0 : DO j = 1, SIZE(gci%lg3x3)
197 : ! Reset langrange multiplier
198 0 : gci%lg3x3(j)%lambda = 0.0_dp
199 : END DO
200 : END IF
201 0 : IF (ASSOCIATED(gci%lg4x6)) THEN
202 0 : DO j = 1, SIZE(gci%lg4x6)
203 : ! Reset langrange multiplier
204 0 : gci%lg4x6(j)%lambda = 0.0_dp
205 : END DO
206 : END IF
207 : END IF
208 : END IF
209 : END IF
210 :
211 24 : IF (has_pos) THEN
212 24 : CALL update_particle_set(particles%els, force_env%para_env, pos=my_pos)
213 24 : DEALLOCATE (my_pos)
214 : END IF
215 24 : IF (has_vel) THEN
216 24 : CALL update_particle_set(particles%els, force_env%para_env, vel=my_vel)
217 24 : DEALLOCATE (my_vel)
218 : END IF
219 24 : CALL timestop(handle)
220 24 : END SUBROUTINE force_env_shake
221 :
222 : ! **************************************************************************************************
223 : !> \brief perform rattle (enforcing of constraints on velocities)
224 : !> This routine can be easily adapted to performe rattle on whatever
225 : !> other vector different from forces..
226 : !> \param force_env the force env to shake
227 : !> \param dt the dt for shake (if you are not interested in the velocities
228 : !> it can be any positive number)
229 : !> \param shake_tol the tolerance for shake
230 : !> \param log_unit if >0 then some information on the shake is printed,
231 : !> defaults to -1
232 : !> \param lagrange_mult ...
233 : !> \param dump_lm ...
234 : !> \param vel ...
235 : !> \param reset ...
236 : !> \author tlaino
237 : ! **************************************************************************************************
238 48 : SUBROUTINE force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
239 24 : vel, reset)
240 :
241 : TYPE(force_env_type), POINTER :: force_env
242 : REAL(kind=dp), INTENT(in), OPTIONAL :: dt
243 : REAL(kind=dp), INTENT(in) :: shake_tol
244 : INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult
245 : LOGICAL, INTENT(IN), OPTIONAL :: dump_lm
246 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
247 : OPTIONAL, TARGET :: vel
248 : LOGICAL, INTENT(IN), OPTIONAL :: reset
249 :
250 : CHARACTER(len=*), PARAMETER :: routineN = 'force_env_rattle'
251 :
252 : INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, &
253 : my_log_unit, nparticle_kind, nparticle_local
254 : LOGICAL :: has_vel, my_dump_lm
255 : REAL(KIND=dp) :: mydt
256 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_vel
257 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
258 : TYPE(cell_type), POINTER :: cell
259 : TYPE(cp_subsys_type), POINTER :: subsys
260 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
261 : TYPE(global_constraint_type), POINTER :: gci
262 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
263 : TYPE(molecule_list_type), POINTER :: molecules
264 : TYPE(particle_list_type), POINTER :: particles
265 :
266 24 : CALL timeset(routineN, handle)
267 24 : CPASSERT(ASSOCIATED(force_env))
268 24 : CPASSERT(force_env%ref_count > 0)
269 24 : my_log_unit = -1
270 24 : IF (PRESENT(log_unit)) my_log_unit = log_unit
271 24 : my_lagrange_mult = -1
272 24 : IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
273 24 : my_dump_lm = .FALSE.
274 24 : IF (PRESENT(dump_lm)) my_dump_lm = dump_lm
275 24 : NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
276 24 : my_vel)
277 24 : IF (PRESENT(vel)) my_vel => vel
278 24 : mydt = 0.1_dp
279 24 : IF (PRESENT(dt)) mydt = dt
280 24 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
281 : CALL cp_subsys_get(subsys, &
282 : atomic_kinds=atomic_kinds, &
283 : local_molecules=local_molecules, &
284 : local_particles=local_particles, &
285 : molecules=molecules, &
286 : molecule_kinds=molecule_kinds, &
287 : particles=particles, &
288 24 : gci=gci)
289 24 : nparticle_kind = atomic_kinds%n_els
290 24 : has_vel = .FALSE.
291 24 : IF (.NOT. ASSOCIATED(my_vel)) THEN
292 24 : has_vel = .TRUE.
293 72 : ALLOCATE (my_vel(3, particles%n_els))
294 7504 : my_vel = 0.0_dp
295 112 : DO iparticle_kind = 1, nparticle_kind
296 88 : nparticle_local = local_particles%n_el(iparticle_kind)
297 1047 : DO iparticle_local = 1, nparticle_local
298 935 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
299 3828 : my_vel(:, iparticle) = particles%els(iparticle)%v(:)
300 : END DO
301 : END DO
302 : END IF
303 :
304 : CALL rattle_control(gci=gci, local_molecules=local_molecules, &
305 : molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, &
306 : particle_set=particles%els, vel=my_vel, dt=mydt, &
307 : rattle_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, &
308 : dump_lm=my_dump_lm, cell=cell, group=force_env%para_env, &
309 24 : local_particles=local_particles)
310 :
311 : ! Possibly reset the lagrange multipliers
312 24 : IF (PRESENT(reset)) THEN
313 24 : IF (reset) THEN
314 : ! Reset Intramolecular constraints
315 500 : DO i = 1, SIZE(molecules%els)
316 476 : IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN
317 28 : DO j = 1, SIZE(molecules%els(i)%lci%lcolv)
318 : ! Reset langrange multiplier
319 28 : molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
320 : END DO
321 : END IF
322 476 : IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN
323 890 : DO j = 1, SIZE(molecules%els(i)%lci%lg3x3)
324 : ! Reset langrange multiplier
325 2216 : molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
326 : END DO
327 : END IF
328 500 : IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN
329 12 : DO j = 1, SIZE(molecules%els(i)%lci%lg4x6)
330 : ! Reset langrange multiplier
331 36 : molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
332 : END DO
333 : END IF
334 : END DO
335 : ! Reset Intermolecular constraints
336 24 : IF (ASSOCIATED(gci)) THEN
337 24 : IF (ASSOCIATED(gci%lcolv)) THEN
338 0 : DO j = 1, SIZE(gci%lcolv)
339 : ! Reset langrange multiplier
340 0 : gci%lcolv(j)%lambda = 0.0_dp
341 : END DO
342 : END IF
343 24 : IF (ASSOCIATED(gci%lg3x3)) THEN
344 0 : DO j = 1, SIZE(gci%lg3x3)
345 : ! Reset langrange multiplier
346 0 : gci%lg3x3(j)%lambda = 0.0_dp
347 : END DO
348 : END IF
349 24 : IF (ASSOCIATED(gci%lg4x6)) THEN
350 0 : DO j = 1, SIZE(gci%lg4x6)
351 : ! Reset langrange multiplier
352 0 : gci%lg4x6(j)%lambda = 0.0_dp
353 : END DO
354 : END IF
355 : END IF
356 : END IF
357 : END IF
358 :
359 24 : IF (has_vel) THEN
360 24 : CALL update_particle_set(particles%els, force_env%para_env, vel=my_vel)
361 : END IF
362 24 : DEALLOCATE (my_vel)
363 24 : CALL timestop(handle)
364 24 : END SUBROUTINE force_env_rattle
365 :
366 : ! **************************************************************************************************
367 : !> \brief Rescale forces if requested
368 : !> \param force_env the force env to shake
369 : !> \author tlaino
370 : ! **************************************************************************************************
371 197716 : SUBROUTINE rescale_forces(force_env)
372 : TYPE(force_env_type), POINTER :: force_env
373 :
374 : CHARACTER(len=*), PARAMETER :: routineN = 'rescale_forces'
375 :
376 : INTEGER :: handle, iparticle
377 : LOGICAL :: explicit
378 : REAL(KIND=dp) :: force(3), max_value, mod_force
379 : TYPE(cp_subsys_type), POINTER :: subsys
380 : TYPE(particle_list_type), POINTER :: particles
381 : TYPE(section_vals_type), POINTER :: rescale_force_section
382 :
383 98858 : CALL timeset(routineN, handle)
384 98858 : CPASSERT(ASSOCIATED(force_env))
385 98858 : CPASSERT(force_env%ref_count > 0)
386 98858 : rescale_force_section => section_vals_get_subs_vals(force_env%force_env_section, "RESCALE_FORCES")
387 98858 : CALL section_vals_get(rescale_force_section, explicit=explicit)
388 98858 : IF (explicit) THEN
389 224 : CALL section_vals_val_get(rescale_force_section, "MAX_FORCE", r_val=max_value)
390 224 : CALL force_env_get(force_env, subsys=subsys)
391 224 : CALL cp_subsys_get(subsys, particles=particles)
392 36262 : DO iparticle = 1, SIZE(particles%els)
393 144152 : force = particles%els(iparticle)%f(:)
394 144152 : mod_force = SQRT(DOT_PRODUCT(force, force))
395 36262 : IF ((mod_force > max_value) .AND. (mod_force /= 0.0_dp)) THEN
396 10208 : force = force/mod_force*max_value
397 10208 : particles%els(iparticle)%f(:) = force
398 : END IF
399 : END DO
400 : END IF
401 98858 : CALL timestop(handle)
402 98858 : END SUBROUTINE rescale_forces
403 :
404 : ! **************************************************************************************************
405 : !> \brief Write forces
406 : !>
407 : !> \param particles ...
408 : !> \param iw ...
409 : !> \param label ...
410 : !> \param ndigits ...
411 : !> \param unit_string ...
412 : !> \param total_force ...
413 : !> \param grand_total_force ...
414 : !> \param zero_force_core_shell_atom ...
415 : !> \author MK (06.09.2010)
416 : ! **************************************************************************************************
417 1820 : SUBROUTINE write_forces(particles, iw, label, ndigits, unit_string, total_force, &
418 : grand_total_force, zero_force_core_shell_atom)
419 :
420 : TYPE(particle_list_type), POINTER :: particles
421 : INTEGER, INTENT(IN) :: iw
422 : CHARACTER(LEN=*), INTENT(IN) :: label
423 : INTEGER, INTENT(IN) :: ndigits
424 : CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
425 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: total_force
426 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
427 : OPTIONAL :: grand_total_force
428 : LOGICAL, INTENT(IN), OPTIONAL :: zero_force_core_shell_atom
429 :
430 : CHARACTER(LEN=18) :: fmtstr4
431 : CHARACTER(LEN=29) :: fmtstr3
432 : CHARACTER(LEN=38) :: fmtstr2
433 : CHARACTER(LEN=49) :: fmtstr1
434 : CHARACTER(LEN=7) :: tag
435 1820 : CHARACTER(LEN=LEN_TRIM(label)) :: lc_label
436 : INTEGER :: i, iparticle, n
437 : LOGICAL :: zero_force
438 : REAL(KIND=dp) :: fconv
439 : REAL(KIND=dp), DIMENSION(3) :: f
440 :
441 1820 : IF (iw > 0) THEN
442 1820 : CPASSERT(ASSOCIATED(particles))
443 1820 : tag = "FORCES|"
444 1820 : lc_label = TRIM(label)
445 1820 : CALL lowercase(lc_label)
446 1820 : n = MIN(MAX(1, ndigits), 20)
447 1820 : fmtstr1 = "(/,T2,A,1X,A,/,T2,A,3X,A,T20,A3,2( X,A3), X,A3)"
448 1820 : WRITE (UNIT=fmtstr1(35:36), FMT="(I2)") n + 5
449 1820 : WRITE (UNIT=fmtstr1(43:44), FMT="(I2)") n + 6
450 1820 : fmtstr2 = "(T2,A,I7,T16,3(1X,ES . ),2X,ES . )"
451 1820 : WRITE (UNIT=fmtstr2(21:22), FMT="(I2)") n + 7
452 1820 : WRITE (UNIT=fmtstr2(24:25), FMT="(I2)") n
453 1820 : WRITE (UNIT=fmtstr2(33:34), FMT="(I2)") n + 7
454 1820 : WRITE (UNIT=fmtstr2(36:37), FMT="(I2)") n
455 1820 : fmtstr3 = "(T2,A,T16,3(1X,ES . ))"
456 1820 : WRITE (UNIT=fmtstr3(18:19), FMT="(I2)") n + 7
457 1820 : WRITE (UNIT=fmtstr3(21:22), FMT="(I2)") n
458 1820 : fmtstr4 = "(T2,A,T ,ES . )"
459 1820 : WRITE (UNIT=fmtstr4(8:9), FMT="(I2)") 3*(n + 8) + 18
460 1820 : WRITE (UNIT=fmtstr4(13:14), FMT="(I2)") n + 7
461 1820 : WRITE (UNIT=fmtstr4(16:17), FMT="(I2)") n
462 1820 : IF (PRESENT(zero_force_core_shell_atom)) THEN
463 501 : zero_force = zero_force_core_shell_atom
464 : ELSE
465 : zero_force = .FALSE.
466 : END IF
467 1820 : fconv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_string))
468 : WRITE (UNIT=iw, FMT=fmtstr1) &
469 1820 : tag, label//" forces ["//TRIM(ADJUSTL(unit_string))//"]", &
470 3640 : tag, "Atom", " x ", " y ", " z ", "|f|"
471 1820 : total_force(1:3) = 0.0_dp
472 188391 : DO iparticle = 1, particles%n_els
473 186571 : IF (particles%els(iparticle)%atom_index /= 0) THEN
474 186571 : i = particles%els(iparticle)%atom_index
475 : ELSE
476 0 : i = iparticle
477 : END IF
478 186571 : IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0)) THEN
479 49176 : f(1:3) = 0.0_dp
480 : ELSE
481 549580 : f(1:3) = particles%els(iparticle)%f(1:3)*fconv
482 : END IF
483 : WRITE (UNIT=iw, FMT=fmtstr2) &
484 932855 : tag, i, f(1:3), SQRT(SUM(f(1:3)**2))
485 748104 : total_force(1:3) = total_force(1:3) + f(1:3)
486 : END DO
487 : WRITE (UNIT=iw, FMT=fmtstr3) &
488 1820 : tag//" Sum", total_force(1:3)
489 : WRITE (UNIT=iw, FMT=fmtstr4) &
490 1820 : tag//" Total "//TRIM(ADJUSTL(lc_label))//" force", &
491 9100 : SQRT(SUM(total_force(1:3)**2))
492 : END IF
493 :
494 1820 : IF (PRESENT(grand_total_force)) THEN
495 668 : grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
496 167 : WRITE (UNIT=iw, FMT="(A)") ""
497 : WRITE (UNIT=iw, FMT=fmtstr4) &
498 167 : tag//" Grand total force ["//TRIM(ADJUSTL(unit_string))//"]", &
499 835 : SQRT(SUM(grand_total_force(1:3)**2))
500 : END IF
501 :
502 1820 : END SUBROUTINE write_forces
503 :
504 : ! **************************************************************************************************
505 : !> \brief Write the atomic coordinates, types, and energies
506 : !> \param iounit ...
507 : !> \param particles ...
508 : !> \param atener ...
509 : !> \param label ...
510 : !> \date 05.06.2023
511 : !> \author JGH
512 : !> \version 1.0
513 : ! **************************************************************************************************
514 489 : SUBROUTINE write_atener(iounit, particles, atener, label)
515 :
516 : INTEGER, INTENT(IN) :: iounit
517 : TYPE(particle_list_type), POINTER :: particles
518 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: atener
519 : CHARACTER(LEN=*), INTENT(IN) :: label
520 :
521 : INTEGER :: i, natom
522 :
523 489 : IF (iounit > 0) THEN
524 489 : WRITE (UNIT=iounit, FMT="(/,T2,A)") TRIM(label)
525 : WRITE (UNIT=iounit, FMT="(T4,A,T30,A,T42,A,T54,A,T69,A)") &
526 489 : "Atom Element", "X", "Y", "Z", "Energy[a.u.]"
527 489 : natom = particles%n_els
528 17788 : DO i = 1, natom
529 17299 : WRITE (iounit, "(I6,T12,A2,T24,3F12.6,F21.10)") i, &
530 17299 : TRIM(ADJUSTL(particles%els(i)%atomic_kind%element_symbol)), &
531 86984 : particles%els(i)%r(1:3)*angstrom, atener(i)
532 : END DO
533 489 : WRITE (iounit, "(A)") ""
534 : END IF
535 :
536 489 : END SUBROUTINE write_atener
537 :
538 : END MODULE force_env_utils
|
