(git:97501a3)
Loading...
Searching...
No Matches
motion_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 Output Utilities for MOTION_SECTION
10!> \author Teodoro Laino [tlaino] - University of Zurich
11!> \date 02.2008
12! **************************************************************************************************
13MODULE motion_utils
14
15 USE cell_types, ONLY: cell_type
16 USE cp2k_info, ONLY: compile_revision,&
17 cp2k_version,&
18 r_host_name,&
19 r_user_name
20 USE cp_log_handling, ONLY: cp_get_default_logger,&
21 cp_logger_type
22 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
23 cp_print_key_unit_nr
24 USE cp_subsys_types, ONLY: cp_subsys_get,&
25 cp_subsys_type
26 USE cp_units, ONLY: cp_unit_from_cp2k
27 USE force_env_types, ONLY: force_env_get,&
28 force_env_type
29 USE input_constants, ONLY: dump_atomic,&
30 dump_dcd,&
31 dump_dcd_aligned_cell,&
32 dump_pdb,&
33 dump_xmol
34 USE input_section_types, ONLY: section_get_ival,&
35 section_vals_get,&
36 section_vals_get_subs_vals,&
37 section_vals_type,&
38 section_vals_val_get
39 USE kinds, ONLY: default_string_length,&
40 dp,&
41 sp
42 USE machine, ONLY: m_flush,&
43 m_timestamp,&
44 timestamp_length
45 USE mathlib, ONLY: diamat_all
46 USE particle_list_types, ONLY: particle_list_type
47 USE particle_methods, ONLY: write_particle_coordinates
48 USE particle_types, ONLY: particle_type
49 USE physcon, ONLY: angstrom
50 USE virial_types, ONLY: virial_type
51#include "./base/base_uses.f90"
52
53 IMPLICIT NONE
54
55 PRIVATE
56
57 PUBLIC :: write_trajectory, write_stress_tensor_to_file, write_simulation_cell, &
58 get_output_format, rot_ana
59
60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'motion_utils'
61 REAL(kind=dp), PARAMETER, PUBLIC :: thrs_motion = 5.0e-10_dp
62
63CONTAINS
64
65! **************************************************************************************************
66!> \brief Performs an analysis of the principal inertia axis
67!> Getting back the generators of the translating and
68!> rotating frame
69!> \param particles ...
70!> \param mat ...
71!> \param dof ...
72!> \param print_section ...
73!> \param keep_rotations ...
74!> \param mass_weighted ...
75!> \param natoms ...
76!> \param rot_dof ...
77!> \param inertia ...
78!> \author Teodoro Laino 08.2006
79! **************************************************************************************************
80 SUBROUTINE rot_ana(particles, mat, dof, print_section, keep_rotations, mass_weighted, &
81 natoms, rot_dof, inertia)
82 TYPE(particle_type), DIMENSION(:), POINTER :: particles
83 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: mat
84 INTEGER, INTENT(OUT) :: dof
85 TYPE(section_vals_type), POINTER :: print_section
86 LOGICAL, INTENT(IN) :: keep_rotations, mass_weighted
87 INTEGER, INTENT(IN) :: natoms
88 INTEGER, INTENT(OUT), OPTIONAL :: rot_dof
89 REAL(kind=dp), INTENT(OUT), OPTIONAL :: inertia(3)
90
91 CHARACTER(len=*), PARAMETER :: routinen = 'rot_ana'
92
93 INTEGER :: handle, i, iparticle, iseq, iw, j, k, &
94 lrot(3)
95 LOGICAL :: present_mat
96 REAL(kind=dp) :: cp(3), ip(3, 3), ip_eigval(3), mass, &
97 masst, norm, rcom(3), rm(3)
98 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rot, tr
99 TYPE(cp_logger_type), POINTER :: logger
100
101 CALL timeset(routinen, handle)
102 logger => cp_get_default_logger()
103 present_mat = PRESENT(mat)
104 cpassert(ASSOCIATED(particles))
105 IF (present_mat) THEN
106 cpassert(.NOT. ASSOCIATED(mat))
107 END IF
108 IF (.NOT. keep_rotations) THEN
109 rcom = 0.0_dp
110 masst = 0.0_dp
111 ! Center of mass
112 DO iparticle = 1, natoms
113 mass = 1.0_dp
114 IF (mass_weighted) mass = particles(iparticle)%atomic_kind%mass
115 cpassert(mass >= 0.0_dp)
116 masst = masst + mass
117 rcom = particles(iparticle)%r*mass + rcom
118 END DO
119 cpassert(masst > 0.0_dp)
120 rcom = rcom/masst
121 ! Intertia Tensor
122 ip = 0.0_dp
123 DO iparticle = 1, natoms
124 mass = 1.0_dp
125 IF (mass_weighted) mass = particles(iparticle)%atomic_kind%mass
126 rm = particles(iparticle)%r - rcom
127 ip(1, 1) = ip(1, 1) + mass*(rm(2)**2 + rm(3)**2)
128 ip(2, 2) = ip(2, 2) + mass*(rm(1)**2 + rm(3)**2)
129 ip(3, 3) = ip(3, 3) + mass*(rm(1)**2 + rm(2)**2)
130 ip(1, 2) = ip(1, 2) - mass*(rm(1)*rm(2))
131 ip(1, 3) = ip(1, 3) - mass*(rm(1)*rm(3))
132 ip(2, 3) = ip(2, 3) - mass*(rm(2)*rm(3))
133 END DO
134 ! Diagonalize the Inertia Tensor
135 CALL diamat_all(ip, ip_eigval)
136 IF (PRESENT(inertia)) inertia = ip_eigval
137 iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
138 IF (iw > 0) THEN
139 WRITE (unit=iw, fmt='(/,T2,A)') &
140 'ROT| Rotational analysis information'
141 WRITE (unit=iw, fmt='(T2,A)') &
142 'ROT| Principal axes and moments of inertia [a.u.]'
143 WRITE (unit=iw, fmt='(T2,A,T14,3(1X,I19))') &
144 'ROT|', 1, 2, 3
145 WRITE (unit=iw, fmt='(T2,A,T21,3(1X,ES19.11))') &
146 'ROT| Eigenvalues', ip_eigval(1:3)
147 WRITE (unit=iw, fmt='(T2,A,T21,3(1X,F19.12))') &
148 'ROT| x', ip(1, 1:3)
149 WRITE (unit=iw, fmt='(T2,A,T21,3(1X,F19.12))') &
150 'ROT| y', ip(2, 1:3)
151 WRITE (unit=iw, fmt='(T2,A,T21,3(1X,F19.12))') &
152 'ROT| z', ip(3, 1:3)
153 END IF
154 CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
155 iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO/COORDINATES", extension=".vibLog")
156 IF (iw > 0) THEN
157 WRITE (unit=iw, fmt='(/,T2,A)') 'ROT| Standard molecule orientation in Angstrom'
158 DO iparticle = 1, natoms
159 WRITE (unit=iw, fmt='(T2,"ROT|",T20,A,T27,3(3X,F15.9))') &
160 trim(particles(iparticle)%atomic_kind%name), &
161 matmul(particles(iparticle)%r, ip)*angstrom
162 END DO
163 END IF
164 CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO/COORDINATES")
165 END IF
166 ! Build up the Translational vectors
167 ALLOCATE (tr(natoms*3, 3))
168 tr = 0.0_dp
169 DO k = 1, 3
170 iseq = 0
171 DO iparticle = 1, natoms
172 mass = 1.0_dp
173 IF (mass_weighted) mass = sqrt(particles(iparticle)%atomic_kind%mass)
174 DO j = 1, 3
175 iseq = iseq + 1
176 IF (j == k) tr(iseq, k) = mass
177 END DO
178 END DO
179 END DO
180 ! Normalize Translations
181 DO i = 1, 3
182 norm = sqrt(dot_product(tr(:, i), tr(:, i)))
183 tr(:, i) = tr(:, i)/norm
184 END DO
185 dof = 3
186 ! Build up the Rotational vectors
187 ALLOCATE (rot(natoms*3, 3))
188 lrot = 0
189 IF (.NOT. keep_rotations) THEN
190 DO iparticle = 1, natoms
191 mass = 1.0_dp
192 IF (mass_weighted) mass = sqrt(particles(iparticle)%atomic_kind%mass)
193 rm = particles(iparticle)%r - rcom
194 cp(1) = rm(1)*ip(1, 1) + rm(2)*ip(2, 1) + rm(3)*ip(3, 1)
195 cp(2) = rm(1)*ip(1, 2) + rm(2)*ip(2, 2) + rm(3)*ip(3, 2)
196 cp(3) = rm(1)*ip(1, 3) + rm(2)*ip(2, 3) + rm(3)*ip(3, 3)
197 ! X Rot
198 rot((iparticle - 1)*3 + 1, 1) = (cp(2)*ip(1, 3) - ip(1, 2)*cp(3))*mass
199 rot((iparticle - 1)*3 + 2, 1) = (cp(2)*ip(2, 3) - ip(2, 2)*cp(3))*mass
200 rot((iparticle - 1)*3 + 3, 1) = (cp(2)*ip(3, 3) - ip(3, 2)*cp(3))*mass
201 ! Y Rot
202 rot((iparticle - 1)*3 + 1, 2) = (cp(3)*ip(1, 1) - ip(1, 3)*cp(1))*mass
203 rot((iparticle - 1)*3 + 2, 2) = (cp(3)*ip(2, 1) - ip(2, 3)*cp(1))*mass
204 rot((iparticle - 1)*3 + 3, 2) = (cp(3)*ip(3, 1) - ip(3, 3)*cp(1))*mass
205 ! Z Rot
206 rot((iparticle - 1)*3 + 1, 3) = (cp(1)*ip(1, 2) - ip(1, 1)*cp(2))*mass
207 rot((iparticle - 1)*3 + 2, 3) = (cp(1)*ip(2, 2) - ip(2, 1)*cp(2))*mass
208 rot((iparticle - 1)*3 + 3, 3) = (cp(1)*ip(3, 2) - ip(3, 1)*cp(2))*mass
209 END DO
210
211 ! Normalize Rotations and count the number of degree of freedom
212 lrot = 1
213 DO i = 1, 3
214 norm = dot_product(rot(:, i), rot(:, i))
215 IF (norm <= thrs_motion) THEN
216 lrot(i) = 0
217 cycle
218 END IF
219 rot(:, i) = rot(:, i)/sqrt(norm)
220 ! Clean Rotational modes for spurious/numerical contamination
221 IF (i < 3) THEN
222 DO j = 1, i
223 rot(:, i + 1) = rot(:, i + 1) - dot_product(rot(:, i + 1), rot(:, j))*rot(:, j)
224 END DO
225 END IF
226 END DO
227 END IF
228 IF (PRESENT(rot_dof)) rot_dof = count(lrot == 1)
229 dof = dof + count(lrot == 1)
230 iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
231 IF (iw > 0) THEN
232 WRITE (iw, '(T2,A,T71,I10)') 'ROT| Number of rotovibrational vectors', dof
233 IF (dof == 5) THEN
234 WRITE (iw, '(T2,A)') &
235 'ROT| Linear molecule detected'
236 END IF
237 IF ((dof == 3) .AND. (.NOT. keep_rotations)) THEN
238 WRITE (iw, '(T2,A)') &
239 'ROT| Single atom detected'
240 END IF
241 END IF
242 CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
243 IF (present_mat) THEN
244 ! Give back the vectors generating the rototranslating Frame
245 ALLOCATE (mat(natoms*3, dof))
246 iseq = 0
247 DO i = 1, 3
248 mat(:, i) = tr(:, i)
249 IF (lrot(i) == 1) THEN
250 iseq = iseq + 1
251 mat(:, 3 + iseq) = rot(:, i)
252 END IF
253 END DO
254 END IF
255 DEALLOCATE (tr)
256 DEALLOCATE (rot)
257 CALL timestop(handle)
258
259 END SUBROUTINE rot_ana
260
261! **************************************************************************************************
262!> \brief Prints the information controlled by the TRAJECTORY section
263!> \param force_env ...
264!> \param root_section ...
265!> \param it ...
266!> \param time ...
267!> \param dtime ...
268!> \param etot ...
269!> \param pk_name ...
270!> \param pos ...
271!> \param act ...
272!> \param middle_name ...
273!> \param particles ...
274!> \param extended_xmol_title ...
275!> \date 02.2008
276!> \author Teodoro Laino [tlaino] - University of Zurich
277!> \version 1.0
278! **************************************************************************************************
279 SUBROUTINE write_trajectory(force_env, root_section, it, time, dtime, etot, pk_name, &
280 pos, act, middle_name, particles, extended_xmol_title)
281 TYPE(force_env_type), POINTER :: force_env
282 TYPE(section_vals_type), POINTER :: root_section
283 INTEGER, INTENT(IN) :: it
284 REAL(kind=dp), INTENT(IN) :: time, dtime, etot
285 CHARACTER(LEN=*), OPTIONAL :: pk_name
286 CHARACTER(LEN=default_string_length), OPTIONAL :: pos, act
287 CHARACTER(LEN=*), OPTIONAL :: middle_name
288 TYPE(particle_list_type), OPTIONAL, POINTER :: particles
289 LOGICAL, INTENT(IN), OPTIONAL :: extended_xmol_title
290
291 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_trajectory'
292
293 CHARACTER(LEN=4) :: id_dcd
294 CHARACTER(LEN=80), DIMENSION(2) :: remark
295 CHARACTER(LEN=default_string_length) :: id_label, id_wpc, my_act, my_ext, &
296 my_form, my_middle, my_pk_name, &
297 my_pos, section_ref, title, unit_str
298 CHARACTER(LEN=timestamp_length) :: timestamp
299 INTEGER :: handle, i, ii, iskip, nat, outformat, &
300 traj_unit
301 INTEGER, POINTER :: force_mixing_indices(:), &
302 force_mixing_labels(:)
303 LOGICAL :: charge_beta, charge_extended, &
304 charge_occup, explicit, &
305 my_extended_xmol_title, new_file, &
306 print_kind
307 REAL(dp), ALLOCATABLE :: fml_array(:)
308 REAL(kind=dp) :: unit_conv
309 TYPE(cell_type), POINTER :: cell
310 TYPE(cp_logger_type), POINTER :: logger
311 TYPE(cp_subsys_type), POINTER :: subsys
312 TYPE(particle_list_type), POINTER :: my_particles
313 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
314 TYPE(section_vals_type), POINTER :: force_env_section, &
315 force_mixing_restart_section
316
317 CALL timeset(routinen, handle)
318
319 NULLIFY (logger, cell, subsys, my_particles, particle_set)
320 logger => cp_get_default_logger()
321 id_label = logger%iter_info%level_name(logger%iter_info%n_rlevel)
322 my_pos = "APPEND"
323 my_act = "WRITE"
324 my_middle = "pos"
325 my_pk_name = "TRAJECTORY"
326 IF (PRESENT(middle_name)) my_middle = middle_name
327 IF (PRESENT(pos)) my_pos = pos
328 IF (PRESENT(act)) my_act = act
329 IF (PRESENT(pk_name)) my_pk_name = pk_name
330
331 SELECT CASE (trim(my_pk_name))
332 CASE ("TRAJECTORY", "SHELL_TRAJECTORY", "CORE_TRAJECTORY")
333 id_dcd = "CORD"
334 id_wpc = "POS"
335 CASE ("VELOCITIES", "SHELL_VELOCITIES", "CORE_VELOCITIES")
336 id_dcd = "VEL "
337 id_wpc = "VEL"
338 CASE ("FORCES", "SHELL_FORCES", "CORE_FORCES")
339 id_dcd = "FRC "
340 id_wpc = "FORCE"
341 CASE ("FORCE_MIXING_LABELS")
342 id_dcd = "FML "
343 id_wpc = "FORCE_MIXING_LABELS"
344 CASE DEFAULT
345 cpabort("")
346 END SELECT
347
348 charge_occup = .false.
349 charge_beta = .false.
350 charge_extended = .false.
351 print_kind = .false.
352
353 CALL force_env_get(force_env, cell=cell, subsys=subsys)
354 IF (PRESENT(particles)) THEN
355 cpassert(ASSOCIATED(particles))
356 my_particles => particles
357 ELSE
358 CALL cp_subsys_get(subsys=subsys, particles=my_particles)
359 END IF
360 particle_set => my_particles%els
361 nat = my_particles%n_els
362
363 ! Gather units of measure for output (if available)
364 IF (trim(my_pk_name) /= "FORCE_MIXING_LABELS") THEN
365 CALL section_vals_val_get(root_section, "MOTION%PRINT%"//trim(my_pk_name)//"%UNIT", &
366 c_val=unit_str)
367 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
368 END IF
369
370 ! Get the output format
371 CALL get_output_format(root_section, "MOTION%PRINT%"//trim(my_pk_name), my_form, my_ext)
372 traj_unit = cp_print_key_unit_nr(logger, root_section, "MOTION%PRINT%"//trim(my_pk_name), &
373 extension=my_ext, file_position=my_pos, file_action=my_act, &
374 file_form=my_form, middle_name=trim(my_middle), is_new_file=new_file)
375 IF (traj_unit > 0) THEN
376 CALL section_vals_val_get(root_section, "MOTION%PRINT%"//trim(my_pk_name)//"%FORMAT", &
377 i_val=outformat)
378 title = ""
379 SELECT CASE (outformat)
380 CASE (dump_dcd, dump_dcd_aligned_cell)
381 IF (new_file) THEN
382 !Lets write the header for the coordinate dcd
383 section_ref = "MOTION%PRINT%"//trim(my_pk_name)//"%EACH%"//trim(id_label)
384 iskip = section_get_ival(root_section, trim(section_ref))
385 i = index(cp2k_version, "(") - 1
386 IF (i == -1) i = len(cp2k_version)
387 CALL m_timestamp(timestamp)
388 WRITE (unit=traj_unit) id_dcd, 0, it, iskip, 0, 0, 0, 0, 0, 0, real(dtime, kind=sp), &
389 1, 0, 0, 0, 0, 0, 0, 0, 0, 24
390 remark(1) = "REMARK "//id_dcd//" DCD file created by "//trim(cp2k_version(:i))// &
391 " (revision "//trim(compile_revision)//")"
392 remark(2) = "REMARK "//trim(r_user_name)//"@"//trim(r_host_name)//" "//timestamp(:19)
393 WRITE (unit=traj_unit) SIZE(remark), remark(:)
394 WRITE (unit=traj_unit) nat
395 CALL m_flush(traj_unit)
396 END IF
397 CASE (dump_xmol)
398 my_extended_xmol_title = .false.
399 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
400 l_val=print_kind)
401 IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
402 ! This information can be digested by Molden
403 IF (my_extended_xmol_title) THEN
404 WRITE (unit=title, fmt="(A,I8,A,F12.3,A,F20.10)") &
405 " i = ", it, ", time = ", time, ", E = ", etot
406 ELSE
407 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i = ", it, ", E = ", etot
408 END IF
409 CASE (dump_atomic)
410 ! Do nothing
411 CASE (dump_pdb)
412 IF (id_wpc == "POS") THEN
413 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_OCCUP", &
414 l_val=charge_occup)
415 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_BETA", &
416 l_val=charge_beta)
417 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_EXTENDED", &
418 l_val=charge_extended)
419 i = count((/charge_occup, charge_beta, charge_extended/))
420 IF (i > 1) &
421 cpabort("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected, ")
422 END IF
423 IF (new_file) THEN
424 CALL m_timestamp(timestamp)
425 ! COLUMNS DATA TYPE FIELD DEFINITION
426 ! 1 - 6 Record name "TITLE "
427 ! 9 - 10 Continuation continuation Allows concatenation
428 ! 11 - 70 String title Title of the experiment
429 WRITE (unit=traj_unit, fmt="(A6,T11,A)") &
430 "TITLE ", "PDB file created by "//trim(cp2k_version)//" (revision "//trim(compile_revision)//")", &
431 "AUTHOR", trim(r_user_name)//"@"//trim(r_host_name)//" "//timestamp(:19)
432 END IF
433 my_extended_xmol_title = .false.
434 IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
435 IF (my_extended_xmol_title) THEN
436 WRITE (unit=title, fmt="(A,I0,A,F0.3,A,F0.10)") &
437 "Step ", it, ", time = ", time, ", E = ", etot
438 ELSE
439 WRITE (unit=title, fmt="(A,I0,A,F0.10)") &
440 "Step ", it, ", E = ", etot
441 END IF
442 CASE DEFAULT
443 cpabort("")
444 END SELECT
445 IF (trim(my_pk_name) == "FORCE_MIXING_LABELS") THEN
446 ALLOCATE (fml_array(3*SIZE(particle_set)))
447 fml_array = 0.0_dp
448 CALL force_env_get(force_env, force_env_section=force_env_section)
449 force_mixing_restart_section => section_vals_get_subs_vals(force_env_section, &
450 "QMMM%FORCE_MIXING%RESTART_INFO", &
451 can_return_null=.true.)
452 IF (ASSOCIATED(force_mixing_restart_section)) THEN
453 CALL section_vals_get(force_mixing_restart_section, explicit=explicit)
454 IF (explicit) THEN
455 CALL section_vals_val_get(force_mixing_restart_section, "INDICES", i_vals=force_mixing_indices)
456 CALL section_vals_val_get(force_mixing_restart_section, "LABELS", i_vals=force_mixing_labels)
457 DO i = 1, SIZE(force_mixing_indices)
458 ii = force_mixing_indices(i)
459 cpassert(ii <= SIZE(particle_set))
460 fml_array((ii - 1)*3 + 1:(ii - 1)*3 + 3) = force_mixing_labels(i)
461 END DO
462 END IF
463 END IF
464 CALL write_particle_coordinates(particle_set, traj_unit, outformat, trim(id_wpc), trim(title), cell, &
465 array=fml_array, print_kind=print_kind)
466 DEALLOCATE (fml_array)
467 ELSE
468 CALL write_particle_coordinates(particle_set, traj_unit, outformat, trim(id_wpc), trim(title), cell, &
469 unit_conv=unit_conv, print_kind=print_kind, &
470 charge_occup=charge_occup, &
471 charge_beta=charge_beta, &
472 charge_extended=charge_extended)
473 END IF
474 END IF
475
476 CALL cp_print_key_finished_output(traj_unit, logger, root_section, "MOTION%PRINT%"//trim(my_pk_name))
477
478 CALL timestop(handle)
479
480 END SUBROUTINE write_trajectory
481
482! **************************************************************************************************
483!> \brief Info on the unit to be opened to dump MD informations
484!> \param section ...
485!> \param path ...
486!> \param my_form ...
487!> \param my_ext ...
488!> \author Teodoro Laino - University of Zurich - 07.2007
489! **************************************************************************************************
490 SUBROUTINE get_output_format(section, path, my_form, my_ext)
491
492 TYPE(section_vals_type), POINTER :: section
493 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: path
494 CHARACTER(LEN=*), INTENT(OUT) :: my_form, my_ext
495
496 INTEGER :: output_format
497
498 IF (PRESENT(path)) THEN
499 CALL section_vals_val_get(section, trim(path)//"%FORMAT", i_val=output_format)
500 ELSE
501 CALL section_vals_val_get(section, "FORMAT", i_val=output_format)
502 END IF
503
504 SELECT CASE (output_format)
505 CASE (dump_dcd, dump_dcd_aligned_cell)
506 my_form = "UNFORMATTED"
507 my_ext = ".dcd"
508 CASE (dump_pdb)
509 my_form = "FORMATTED"
510 my_ext = ".pdb"
511 CASE DEFAULT
512 my_form = "FORMATTED"
513 my_ext = ".xyz"
514 END SELECT
515
516 END SUBROUTINE get_output_format
517
518! **************************************************************************************************
519!> \brief Prints the Stress Tensor
520!> \param virial ...
521!> \param cell ...
522!> \param motion_section ...
523!> \param itimes ...
524!> \param time ...
525!> \param pos ...
526!> \param act ...
527!> \date 02.2008
528!> \author Teodoro Laino [tlaino] - University of Zurich
529!> \version 1.0
530! **************************************************************************************************
531 SUBROUTINE write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
532
533 TYPE(virial_type), POINTER :: virial
534 TYPE(cell_type), POINTER :: cell
535 TYPE(section_vals_type), POINTER :: motion_section
536 INTEGER, INTENT(IN) :: itimes
537 REAL(kind=dp), INTENT(IN) :: time
538 CHARACTER(LEN=default_string_length), INTENT(IN), &
539 OPTIONAL :: pos, act
540
541 CHARACTER(LEN=default_string_length) :: my_act, my_pos
542 INTEGER :: output_unit
543 LOGICAL :: new_file
544 REAL(kind=dp), DIMENSION(3, 3) :: pv_total_bar
545 TYPE(cp_logger_type), POINTER :: logger
546
547 NULLIFY (logger)
548 logger => cp_get_default_logger()
549
550 IF (virial%pv_availability) THEN
551 my_pos = "APPEND"
552 my_act = "WRITE"
553 IF (PRESENT(pos)) my_pos = pos
554 IF (PRESENT(act)) my_act = act
555 output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%STRESS", &
556 extension=".stress", file_position=my_pos, &
557 file_action=my_act, file_form="FORMATTED", &
558 is_new_file=new_file)
559 ELSE
560 output_unit = 0
561 END IF
562
563 IF (output_unit > 0) THEN
564 IF (new_file) THEN
565 WRITE (unit=output_unit, fmt='(A,9(12X,A2," [bar]"),6X,A)') &
566 "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
567 END IF
568 pv_total_bar(1, 1) = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
569 pv_total_bar(1, 2) = cp_unit_from_cp2k(virial%pv_total(1, 2)/cell%deth, "bar")
570 pv_total_bar(1, 3) = cp_unit_from_cp2k(virial%pv_total(1, 3)/cell%deth, "bar")
571 pv_total_bar(2, 1) = cp_unit_from_cp2k(virial%pv_total(2, 1)/cell%deth, "bar")
572 pv_total_bar(2, 2) = cp_unit_from_cp2k(virial%pv_total(2, 2)/cell%deth, "bar")
573 pv_total_bar(2, 3) = cp_unit_from_cp2k(virial%pv_total(2, 3)/cell%deth, "bar")
574 pv_total_bar(3, 1) = cp_unit_from_cp2k(virial%pv_total(3, 1)/cell%deth, "bar")
575 pv_total_bar(3, 2) = cp_unit_from_cp2k(virial%pv_total(3, 2)/cell%deth, "bar")
576 pv_total_bar(3, 3) = cp_unit_from_cp2k(virial%pv_total(3, 3)/cell%deth, "bar")
577 WRITE (unit=output_unit, fmt='(I8,F12.3,9(1X,F19.10))') itimes, time, &
578 pv_total_bar(1, 1), pv_total_bar(1, 2), pv_total_bar(1, 3), &
579 pv_total_bar(2, 1), pv_total_bar(2, 2), pv_total_bar(2, 3), &
580 pv_total_bar(3, 1), pv_total_bar(3, 2), pv_total_bar(3, 3)
581 CALL m_flush(output_unit)
582 END IF
583
584 IF (virial%pv_availability) THEN
585 CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
586 "PRINT%STRESS")
587 END IF
588
589 END SUBROUTINE write_stress_tensor_to_file
590
591! **************************************************************************************************
592!> \brief Prints the Simulation Cell
593!> \param cell ...
594!> \param motion_section ...
595!> \param itimes ...
596!> \param time ...
597!> \param pos ...
598!> \param act ...
599!> \date 02.2008
600!> \author Teodoro Laino [tlaino] - University of Zurich
601!> \version 1.0
602! **************************************************************************************************
603 SUBROUTINE write_simulation_cell(cell, motion_section, itimes, time, pos, act)
604
605 TYPE(cell_type), POINTER :: cell
606 TYPE(section_vals_type), POINTER :: motion_section
607 INTEGER, INTENT(IN) :: itimes
608 REAL(kind=dp), INTENT(IN) :: time
609 CHARACTER(LEN=default_string_length), INTENT(IN), &
610 OPTIONAL :: pos, act
611
612 CHARACTER(LEN=default_string_length) :: my_act, my_pos
613 INTEGER :: output_unit
614 LOGICAL :: new_file
615 TYPE(cp_logger_type), POINTER :: logger
616
617 NULLIFY (logger)
618 logger => cp_get_default_logger()
619
620 my_pos = "APPEND"
621 my_act = "WRITE"
622 IF (PRESENT(pos)) my_pos = pos
623 IF (PRESENT(act)) my_act = act
624
625 output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%CELL", &
626 extension=".cell", file_position=my_pos, &
627 file_action=my_act, file_form="FORMATTED", &
628 is_new_file=new_file)
629
630 IF (output_unit > 0) THEN
631 IF (new_file) THEN
632 WRITE (unit=output_unit, fmt='(A,9(7X,A2," [Angstrom]"),6X,A)') &
633 "# Step Time [fs]", "Ax", "Ay", "Az", "Bx", "By", "Bz", "Cx", "Cy", "Cz", &
634 "Volume [Angstrom^3]"
635 END IF
636 WRITE (unit=output_unit, fmt="(I8,F12.3,9(1X,F19.10),1X,F24.10)") itimes, time, &
637 cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom, &
638 cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom, &
639 cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom, &
640 cell%deth*angstrom*angstrom*angstrom
641 CALL m_flush(output_unit)
642 END IF
643
644 CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
645 "PRINT%CELL")
646
647 END SUBROUTINE write_simulation_cell
648
649END MODULE motion_utils
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp2k_info
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:16
cp2k_info::r_host_name
character(len=default_string_length), public r_host_name
Definition cp2k_info.F:68
cp2k_info::compile_revision
character(len= *), parameter, public compile_revision
Definition cp2k_info.F:39
cp2k_info::cp2k_version
character(len= *), parameter, public cp2k_version
Definition cp2k_info.F:43
cp2k_info::r_user_name
character(len=default_string_length), public r_user_name
Definition cp2k_info.F:68
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_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
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
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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::dump_xmol
integer, parameter, public dump_xmol
Definition input_constants.F:376
input_constants::dump_pdb
integer, parameter, public dump_pdb
Definition input_constants.F:376
input_constants::dump_atomic
integer, parameter, public dump_atomic
Definition input_constants.F:376
input_constants::dump_dcd_aligned_cell
integer, parameter, public dump_dcd_aligned_cell
Definition input_constants.F:376
input_constants::dump_dcd
integer, parameter, public dump_dcd
Definition input_constants.F:376
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_get_ival
integer function, public section_get_ival(section_vals, keyword_name)
...
Definition input_section_types.F:980
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
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::timestamp_length
integer, parameter, public timestamp_length
Definition machine.F:53
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:131
machine::m_timestamp
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:380
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:373
motion_utils
Output Utilities for MOTION_SECTION.
Definition motion_utils.F:13
motion_utils::thrs_motion
real(kind=dp), parameter, public thrs_motion
Definition motion_utils.F:61
motion_utils::get_output_format
subroutine, public get_output_format(section, path, my_form, my_ext)
Info on the unit to be opened to dump MD informations.
Definition motion_utils.F:491
motion_utils::write_simulation_cell
subroutine, public write_simulation_cell(cell, motion_section, itimes, time, pos, act)
Prints the Simulation Cell.
Definition motion_utils.F:604
motion_utils::write_trajectory
subroutine, public write_trajectory(force_env, root_section, it, time, dtime, etot, pk_name, pos, act, middle_name, particles, extended_xmol_title)
Prints the information controlled by the TRAJECTORY section.
Definition motion_utils.F:281
motion_utils::rot_ana
subroutine, public rot_ana(particles, mat, dof, print_section, keep_rotations, mass_weighted, natoms, rot_dof, inertia)
Performs an analysis of the principal inertia axis Getting back the generators of the translating and...
Definition motion_utils.F:82
motion_utils::write_stress_tensor_to_file
subroutine, public write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
Prints the Stress Tensor.
Definition motion_utils.F:532
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::write_particle_coordinates
subroutine, public write_particle_coordinates(particle_set, iunit, output_format, content, title, cell, array, unit_conv, charge_occup, charge_beta, charge_extended, print_kind)
Should be able to write a few formats e.g. xmol, and some binary format (dcd) some format can be used...
Definition particle_methods.F:183
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
virial_types
Definition virial_types.F:13
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
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
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
virial_types::virial_type
Definition virial_types.F:26