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