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