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