(git:71c3ab0)
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 = norm2(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=5) :: pbc_str
297 CHARACTER(LEN=80), DIMENSION(2) :: remark
298 CHARACTER(LEN=default_string_length) :: etot_str, id_extxyz, id_label, id_wpc, my_act, &
299 my_ext, my_form, my_middle, my_pk_name, my_pos, section_ref, step_str, time_str, unit_str
300 CHARACTER(LEN=timestamp_length) :: timestamp
301 INTEGER :: handle, i, ii, iskip, nat, outformat, &
302 traj_unit
303 INTEGER, POINTER :: force_mixing_indices(:), &
304 force_mixing_labels(:)
305 LOGICAL :: charge_beta, charge_extended, &
306 charge_occup, explicit, &
307 my_extended_xmol_title, new_file, &
308 print_kind
309 REAL(dp), ALLOCATABLE :: fml_array(:)
310 REAL(kind=dp) :: unit_conv
311 TYPE(cell_type), POINTER :: cell
312 TYPE(cp_logger_type), POINTER :: logger
313 TYPE(cp_subsys_type), POINTER :: subsys
314 TYPE(particle_list_type), POINTER :: my_particles
315 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
316 TYPE(section_vals_type), POINTER :: force_env_section, &
317 force_mixing_restart_section
318
319 CALL timeset(routinen, handle)
320
321 NULLIFY (logger, cell, subsys, my_particles, particle_set)
322 logger => cp_get_default_logger()
323 id_label = logger%iter_info%level_name(logger%iter_info%n_rlevel)
324 my_pos = "APPEND"
325 my_act = "WRITE"
326 my_middle = "pos"
327 my_pk_name = "TRAJECTORY"
328 IF (PRESENT(middle_name)) my_middle = middle_name
329 IF (PRESENT(pos)) my_pos = pos
330 IF (PRESENT(act)) my_act = act
331 IF (PRESENT(pk_name)) my_pk_name = pk_name
332
333 SELECT CASE (trim(my_pk_name))
334 CASE ("TRAJECTORY", "SHELL_TRAJECTORY", "CORE_TRAJECTORY")
335 id_dcd = "CORD"
336 id_wpc = "POS"
337 id_extxyz = "pos"
338 CASE ("VELOCITIES", "SHELL_VELOCITIES", "CORE_VELOCITIES")
339 id_dcd = "VEL "
340 id_wpc = "VEL"
341 id_extxyz = "velo"
342 CASE ("FORCES", "SHELL_FORCES", "CORE_FORCES")
343 id_dcd = "FRC "
344 id_wpc = "FORCE"
345 id_extxyz = "force"
346 CASE ("FORCE_MIXING_LABELS")
347 id_dcd = "FML "
348 id_wpc = "FORCE_MIXING_LABELS"
349 id_extxyz = "force_mixing_label" ! non-standard
350 CASE DEFAULT
351 CALL cp_abort(__location__, &
352 "<TRAJECTORY>, <VELOCITIES>, <FORCES>, "// &
353 "<FORCE_MIXING_LABELS> are supported as "// &
354 "the <my_pk_name> for write_trajectory, "// &
355 "found unknown option "// &
356 "<"//trim(my_pk_name)//">")
357 END SELECT
358
359 charge_occup = .false.
360 charge_beta = .false.
361 charge_extended = .false.
362 print_kind = .false.
363
364 CALL force_env_get(force_env, cell=cell, subsys=subsys)
365 IF (PRESENT(particles)) THEN
366 cpassert(ASSOCIATED(particles))
367 my_particles => particles
368 ELSE
369 CALL cp_subsys_get(subsys=subsys, particles=my_particles)
370 END IF
371 particle_set => my_particles%els
372 nat = my_particles%n_els
373 pbc_str = "F F F"
374 IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
375 IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
376 IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
377
378 ! Gather units of measure for output (if available)
379 IF (trim(my_pk_name) /= "FORCE_MIXING_LABELS") THEN
380 CALL section_vals_val_get(root_section, "MOTION%PRINT%"//trim(my_pk_name)//"%UNIT", &
381 c_val=unit_str)
382 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
383 END IF
384
385 ! Get the output format
386 CALL get_output_format(root_section, "MOTION%PRINT%"//trim(my_pk_name), my_form, my_ext)
387 traj_unit = cp_print_key_unit_nr(logger, root_section, "MOTION%PRINT%"//trim(my_pk_name), &
388 extension=my_ext, file_position=my_pos, file_action=my_act, &
389 file_form=my_form, middle_name=trim(my_middle), is_new_file=new_file)
390 IF (traj_unit > 0) THEN
391 CALL section_vals_val_get(root_section, "MOTION%PRINT%"//trim(my_pk_name)//"%FORMAT", &
392 i_val=outformat)
393 title = ""
394 SELECT CASE (outformat)
395 CASE (dump_dcd, dump_dcd_aligned_cell)
396 IF (new_file) THEN
397 !Lets write the header for the coordinate dcd
398 section_ref = "MOTION%PRINT%"//trim(my_pk_name)//"%EACH%"//trim(id_label)
399 iskip = section_get_ival(root_section, trim(section_ref))
400 i = index(cp2k_version, "(") - 1
401 IF (i == -1) i = len(cp2k_version)
402 CALL m_timestamp(timestamp)
403 WRITE (unit=traj_unit) id_dcd, 0, it, iskip, 0, 0, 0, 0, 0, 0, real(dtime, kind=sp), &
404 1, 0, 0, 0, 0, 0, 0, 0, 0, 24
405 remark(1) = "REMARK "//id_dcd//" DCD file created by "//trim(cp2k_version(:i))// &
406 " (revision "//trim(compile_revision)//")"
407 remark(2) = "REMARK "//trim(r_user_name)//"@"//trim(r_host_name)//" "//timestamp(:19)
408 WRITE (unit=traj_unit) SIZE(remark), remark(:)
409 WRITE (unit=traj_unit) nat
410 CALL m_flush(traj_unit)
411 END IF
412 CASE (dump_xmol)
413 my_extended_xmol_title = .false.
414 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
415 l_val=print_kind)
416 IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
417 ! This information can be digested by Molden
418 IF (my_extended_xmol_title) THEN
419 WRITE (unit=title, fmt="(A,I8,A,F12.3,A,F20.10)") &
420 " i = ", it, ", time = ", time, ", E = ", etot
421 ELSE
422 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i = ", it, ", E = ", etot
423 END IF
424 CASE (dump_extxyz)
425 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", l_val=print_kind)
426 WRITE (unit=cell_str, fmt="(9(1X,F19.10))") cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom
427 WRITE (unit=step_str, fmt="(I8)") it
428 WRITE (unit=time_str, fmt="(F12.3)") time
429 WRITE (unit=etot_str, fmt="(F20.10)") etot
430 WRITE (unit=title, fmt="(A)") &
431 'Lattice="'//trim(adjustl(cell_str))//'"'// &
432 ' Properties=species:S:1:'//trim(id_extxyz)//':R:3'// &
433 ' pbc="'//pbc_str//'"'// &
434 ' Step='//trim(adjustl(step_str))// &
435 ' Time='//trim(adjustl(time_str))// &
436 ' Energy='//trim(adjustl(etot_str))
437 CASE (dump_atomic)
438 ! Do nothing
439 CASE (dump_pdb)
440 IF (id_wpc == "POS") THEN
441 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_OCCUP", &
442 l_val=charge_occup)
443 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_BETA", &
444 l_val=charge_beta)
445 CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_EXTENDED", &
446 l_val=charge_extended)
447 i = count([charge_occup, charge_beta, charge_extended])
448 IF (i > 1) THEN
449 cpabort("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected, ")
450 END IF
451 END IF
452 IF (new_file) THEN
453 CALL m_timestamp(timestamp)
454 ! COLUMNS DATA TYPE FIELD DEFINITION
455 ! 1 - 6 Record name "TITLE "
456 ! 9 - 10 Continuation continuation Allows concatenation
457 ! 11 - 70 String title Title of the experiment
458 WRITE (unit=traj_unit, fmt="(A6,T11,A)") &
459 "TITLE ", "PDB file created by "//trim(cp2k_version)//" (revision "//trim(compile_revision)//")", &
460 "AUTHOR", trim(r_user_name)//"@"//trim(r_host_name)//" "//timestamp(:19)
461 END IF
462 my_extended_xmol_title = .false.
463 IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
464 IF (my_extended_xmol_title) THEN
465 WRITE (unit=title, fmt="(A,I0,A,F0.3,A,F0.10)") &
466 "Step ", it, ", time = ", time, ", E = ", etot
467 ELSE
468 WRITE (unit=title, fmt="(A,I0,A,F0.10)") &
469 "Step ", it, ", E = ", etot
470 END IF
471 CASE DEFAULT
472 cpabort("Unknown output format")
473 END SELECT
474 IF (trim(my_pk_name) == "FORCE_MIXING_LABELS") THEN
475 ALLOCATE (fml_array(3*SIZE(particle_set)))
476 fml_array = 0.0_dp
477 CALL force_env_get(force_env, force_env_section=force_env_section)
478 force_mixing_restart_section => section_vals_get_subs_vals(force_env_section, &
479 "QMMM%FORCE_MIXING%RESTART_INFO", &
480 can_return_null=.true.)
481 IF (ASSOCIATED(force_mixing_restart_section)) THEN
482 CALL section_vals_get(force_mixing_restart_section, explicit=explicit)
483 IF (explicit) THEN
484 CALL section_vals_val_get(force_mixing_restart_section, "INDICES", i_vals=force_mixing_indices)
485 CALL section_vals_val_get(force_mixing_restart_section, "LABELS", i_vals=force_mixing_labels)
486 DO i = 1, SIZE(force_mixing_indices)
487 ii = force_mixing_indices(i)
488 cpassert(ii <= SIZE(particle_set))
489 fml_array((ii - 1)*3 + 1:(ii - 1)*3 + 3) = force_mixing_labels(i)
490 END DO
491 END IF
492 END IF
493 CALL write_particle_coordinates(particle_set, traj_unit, outformat, trim(id_wpc), trim(title), cell, &
494 array=fml_array, print_kind=print_kind)
495 DEALLOCATE (fml_array)
496 ELSE
497 CALL write_particle_coordinates(particle_set, traj_unit, outformat, trim(id_wpc), trim(title), cell, &
498 unit_conv=unit_conv, print_kind=print_kind, &
499 charge_occup=charge_occup, &
500 charge_beta=charge_beta, &
501 charge_extended=charge_extended)
502 END IF
503 END IF
504
505 CALL cp_print_key_finished_output(traj_unit, logger, root_section, "MOTION%PRINT%"//trim(my_pk_name))
506
507 CALL timestop(handle)
508
509 END SUBROUTINE write_trajectory
510
511! **************************************************************************************************
512!> \brief Info on the unit to be opened to dump MD informations
513!> \param section ...
514!> \param path ...
515!> \param my_form ...
516!> \param my_ext ...
517!> \author Teodoro Laino - University of Zurich - 07.2007
518! **************************************************************************************************
519 SUBROUTINE get_output_format(section, path, my_form, my_ext)
520
521 TYPE(section_vals_type), POINTER :: section
522 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: path
523 CHARACTER(LEN=*), INTENT(OUT) :: my_form, my_ext
524
525 INTEGER :: output_format
526
527 IF (PRESENT(path)) THEN
528 CALL section_vals_val_get(section, trim(path)//"%FORMAT", i_val=output_format)
529 ELSE
530 CALL section_vals_val_get(section, "FORMAT", i_val=output_format)
531 END IF
532
533 SELECT CASE (output_format)
534 CASE (dump_dcd, dump_dcd_aligned_cell)
535 my_form = "UNFORMATTED"
536 my_ext = ".dcd"
537 CASE (dump_pdb)
538 my_form = "FORMATTED"
539 my_ext = ".pdb"
540 CASE DEFAULT
541 my_form = "FORMATTED"
542 my_ext = ".xyz"
543 END SELECT
544
545 END SUBROUTINE get_output_format
546
547! **************************************************************************************************
548!> \brief Prints the Stress Tensor
549!> \param virial ...
550!> \param cell ...
551!> \param motion_section ...
552!> \param itimes ...
553!> \param time ...
554!> \param pos ...
555!> \param act ...
556!> \date 02.2008
557!> \author Teodoro Laino [tlaino] - University of Zurich
558!> \version 1.0
559! **************************************************************************************************
560 SUBROUTINE write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
561
562 TYPE(virial_type), POINTER :: virial
563 TYPE(cell_type), POINTER :: cell
564 TYPE(section_vals_type), POINTER :: motion_section
565 INTEGER, INTENT(IN) :: itimes
566 REAL(kind=dp), INTENT(IN) :: time
567 CHARACTER(LEN=default_string_length), INTENT(IN), &
568 OPTIONAL :: pos, act
569
570 CHARACTER(LEN=default_string_length) :: my_act, my_pos
571 INTEGER :: output_unit
572 LOGICAL :: new_file
573 REAL(kind=dp), DIMENSION(3, 3) :: pv_total_bar
574 TYPE(cp_logger_type), POINTER :: logger
575
576 NULLIFY (logger)
577 logger => cp_get_default_logger()
578
579 IF (virial%pv_availability) THEN
580 my_pos = "APPEND"
581 my_act = "WRITE"
582 IF (PRESENT(pos)) my_pos = pos
583 IF (PRESENT(act)) my_act = act
584 output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%STRESS", &
585 extension=".stress", file_position=my_pos, &
586 file_action=my_act, file_form="FORMATTED", &
587 is_new_file=new_file)
588 ELSE
589 output_unit = 0
590 END IF
591
592 IF (output_unit > 0) THEN
593 IF (new_file) THEN
594 WRITE (unit=output_unit, fmt='(A,9(12X,A2," [bar]"),6X,A)') &
595 "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
596 END IF
597 pv_total_bar(1, 1) = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
598 pv_total_bar(1, 2) = cp_unit_from_cp2k(virial%pv_total(1, 2)/cell%deth, "bar")
599 pv_total_bar(1, 3) = cp_unit_from_cp2k(virial%pv_total(1, 3)/cell%deth, "bar")
600 pv_total_bar(2, 1) = cp_unit_from_cp2k(virial%pv_total(2, 1)/cell%deth, "bar")
601 pv_total_bar(2, 2) = cp_unit_from_cp2k(virial%pv_total(2, 2)/cell%deth, "bar")
602 pv_total_bar(2, 3) = cp_unit_from_cp2k(virial%pv_total(2, 3)/cell%deth, "bar")
603 pv_total_bar(3, 1) = cp_unit_from_cp2k(virial%pv_total(3, 1)/cell%deth, "bar")
604 pv_total_bar(3, 2) = cp_unit_from_cp2k(virial%pv_total(3, 2)/cell%deth, "bar")
605 pv_total_bar(3, 3) = cp_unit_from_cp2k(virial%pv_total(3, 3)/cell%deth, "bar")
606 WRITE (unit=output_unit, fmt='(I8,F12.3,9(1X,F19.10))') itimes, time, &
607 pv_total_bar(1, 1), pv_total_bar(1, 2), pv_total_bar(1, 3), &
608 pv_total_bar(2, 1), pv_total_bar(2, 2), pv_total_bar(2, 3), &
609 pv_total_bar(3, 1), pv_total_bar(3, 2), pv_total_bar(3, 3)
610 CALL m_flush(output_unit)
611 END IF
612
613 IF (virial%pv_availability) THEN
614 CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
615 "PRINT%STRESS")
616 END IF
617
618 END SUBROUTINE write_stress_tensor_to_file
619
620! **************************************************************************************************
621!> \brief Prints the Simulation Cell
622!> \param cell ...
623!> \param motion_section ...
624!> \param itimes ...
625!> \param time ...
626!> \param pos ...
627!> \param act ...
628!> \date 02.2008
629!> \author Teodoro Laino [tlaino] - University of Zurich
630!> \version 1.0
631! **************************************************************************************************
632 SUBROUTINE write_simulation_cell(cell, motion_section, itimes, time, pos, act)
633
634 TYPE(cell_type), POINTER :: cell
635 TYPE(section_vals_type), POINTER :: motion_section
636 INTEGER, INTENT(IN) :: itimes
637 REAL(kind=dp), INTENT(IN) :: time
638 CHARACTER(LEN=default_string_length), INTENT(IN), &
639 OPTIONAL :: pos, act
640
641 CHARACTER(LEN=default_string_length) :: my_act, my_pos
642 INTEGER :: output_unit
643 LOGICAL :: new_file
644 TYPE(cp_logger_type), POINTER :: logger
645
646 NULLIFY (logger)
647 logger => cp_get_default_logger()
648
649 my_pos = "APPEND"
650 my_act = "WRITE"
651 IF (PRESENT(pos)) my_pos = pos
652 IF (PRESENT(act)) my_act = act
653
654 output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%CELL", &
655 extension=".cell", file_position=my_pos, &
656 file_action=my_act, file_form="FORMATTED", &
657 is_new_file=new_file)
658
659 IF (output_unit > 0) THEN
660 IF (new_file) THEN
661 WRITE (unit=output_unit, fmt='(A,9(7X,A2," [Angstrom]"),6X,A)') &
662 "# Step Time [fs]", "Ax", "Ay", "Az", "Bx", "By", "Bz", "Cx", "Cy", "Cz", &
663 "Volume [Angstrom^3]"
664 END IF
665 WRITE (unit=output_unit, fmt="(I8,F12.3,9(1X,F19.10),1X,F24.10)") itimes, time, &
666 cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom, &
667 cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom, &
668 cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom, &
669 cell%deth*angstrom*angstrom*angstrom
670 CALL m_flush(output_unit)
671 END IF
672
673 CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
674 "PRINT%CELL")
675
676 END SUBROUTINE write_simulation_cell
677
678END 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:22
cp2k_info::r_host_name
character(len=default_string_length), public r_host_name
Definition cp2k_info.F:74
cp2k_info::compile_revision
character(len= *), parameter, public compile_revision
Definition cp2k_info.F:45
cp2k_info::cp2k_version
character(len= *), parameter, public cp2k_version
Definition cp2k_info.F:49
cp2k_info::r_user_name
character(len=default_string_length), public r_user_name
Definition cp2k_info.F:74
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:1073
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, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
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:1251
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:344
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::dump_xmol
integer, parameter, public dump_xmol
Definition input_constants.F:386
input_constants::dump_pdb
integer, parameter, public dump_pdb
Definition input_constants.F:386
input_constants::dump_extxyz
integer, parameter, public dump_extxyz
Definition input_constants.F:386
input_constants::dump_atomic
integer, parameter, public dump_atomic
Definition input_constants.F:386
input_constants::dump_dcd_aligned_cell
integer, parameter, public dump_dcd_aligned_cell
Definition input_constants.F:386
input_constants::dump_dcd
integer, parameter, public dump_dcd
Definition input_constants.F:386
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:985
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:735
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:708
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:1052
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:46
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_timestamp
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:381
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:520
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:633
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:561
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:226
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