Line data Source code
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 : ! **************************************************************************************************
13 : MODULE 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 :
64 : CONTAINS
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 2167 : 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 2167 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rot, Tr
100 : TYPE(cp_logger_type), POINTER :: logger
101 :
102 2167 : CALL timeset(routineN, handle)
103 2167 : logger => cp_get_default_logger()
104 2167 : present_mat = PRESENT(mat)
105 2167 : CPASSERT(ASSOCIATED(particles))
106 2167 : IF (present_mat) THEN
107 398 : CPASSERT(.NOT. ASSOCIATED(mat))
108 : END IF
109 2167 : IF (.NOT. keep_rotations) THEN
110 2165 : rcom = 0.0_dp
111 2165 : masst = 0.0_dp
112 : ! Center of mass
113 494583 : DO iparticle = 1, natoms
114 492418 : mass = 1.0_dp
115 492418 : IF (mass_weighted) mass = particles(iparticle)%atomic_kind%mass
116 490592 : CPASSERT(mass >= 0.0_dp)
117 492418 : masst = masst + mass
118 1971837 : rcom = particles(iparticle)%r*mass + rcom
119 : END DO
120 2165 : CPASSERT(masst > 0.0_dp)
121 8660 : rcom = rcom/masst
122 : ! Intertia Tensor
123 2165 : Ip = 0.0_dp
124 494583 : DO iparticle = 1, natoms
125 492418 : mass = 1.0_dp
126 492418 : IF (mass_weighted) mass = particles(iparticle)%atomic_kind%mass
127 1969672 : rm = particles(iparticle)%r - rcom
128 492418 : Ip(1, 1) = Ip(1, 1) + mass*(rm(2)**2 + rm(3)**2)
129 492418 : Ip(2, 2) = Ip(2, 2) + mass*(rm(1)**2 + rm(3)**2)
130 492418 : Ip(3, 3) = Ip(3, 3) + mass*(rm(1)**2 + rm(2)**2)
131 492418 : Ip(1, 2) = Ip(1, 2) - mass*(rm(1)*rm(2))
132 492418 : Ip(1, 3) = Ip(1, 3) - mass*(rm(1)*rm(3))
133 494583 : Ip(2, 3) = Ip(2, 3) - mass*(rm(2)*rm(3))
134 : END DO
135 : ! Diagonalize the Inertia Tensor
136 2165 : CALL diamat_all(Ip, Ip_eigval)
137 2165 : IF (PRESENT(inertia)) inertia = Ip_eigval
138 2165 : iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
139 2165 : IF (iw > 0) THEN
140 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
141 1002 : 'ROT| Rotational analysis information'
142 : WRITE (UNIT=iw, FMT='(T2,A)') &
143 1002 : 'ROT| Principal axes and moments of inertia [a.u.]'
144 : WRITE (UNIT=iw, FMT='(T2,A,T14,3(1X,I19))') &
145 1002 : 'ROT|', 1, 2, 3
146 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
147 1002 : 'ROT| Eigenvalues', Ip_eigval(1:3)
148 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
149 1002 : 'ROT| x', Ip(1, 1:3)
150 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
151 1002 : 'ROT| y', Ip(2, 1:3)
152 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
153 1002 : 'ROT| z', Ip(3, 1:3)
154 : END IF
155 2165 : CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
156 2165 : iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO/COORDINATES", extension=".vibLog")
157 2165 : IF (iw > 0) THEN
158 62 : WRITE (UNIT=iw, FMT='(/,T2,A)') 'ROT| Standard molecule orientation in Angstrom'
159 428 : DO iparticle = 1, natoms
160 : WRITE (UNIT=iw, FMT='(T2,"ROT|",T20,A,T27,3(3X,F15.9))') &
161 366 : TRIM(particles(iparticle)%atomic_kind%name), &
162 6284 : MATMUL(particles(iparticle)%r, Ip)*angstrom
163 : END DO
164 : END IF
165 2165 : CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO/COORDINATES")
166 : END IF
167 : ! Build up the Translational vectors
168 8668 : ALLOCATE (Tr(natoms*3, 3))
169 2167 : Tr = 0.0_dp
170 8668 : DO k = 1, 3
171 : iseq = 0
172 1485940 : DO iparticle = 1, natoms
173 1477272 : mass = 1.0_dp
174 1477272 : IF (mass_weighted) mass = SQRT(particles(iparticle)%atomic_kind%mass)
175 5915589 : DO j = 1, 3
176 4431816 : iseq = iseq + 1
177 5909088 : IF (j == k) Tr(iseq, k) = mass
178 : END DO
179 : END DO
180 : END DO
181 : ! Normalize Translations
182 8668 : DO i = 1, 3
183 4438317 : norm = SQRT(DOT_PRODUCT(Tr(:, i), Tr(:, i)))
184 4440484 : Tr(:, i) = Tr(:, i)/norm
185 : END DO
186 2167 : dof = 3
187 : ! Build up the Rotational vectors
188 4334 : ALLOCATE (Rot(natoms*3, 3))
189 2167 : lrot = 0
190 2167 : IF (.NOT. keep_rotations) THEN
191 494583 : DO iparticle = 1, natoms
192 492418 : mass = 1.0_dp
193 492418 : IF (mass_weighted) mass = SQRT(particles(iparticle)%atomic_kind%mass)
194 1969672 : rm = particles(iparticle)%r - rcom
195 492418 : cp(1) = rm(1)*Ip(1, 1) + rm(2)*Ip(2, 1) + rm(3)*Ip(3, 1)
196 492418 : cp(2) = rm(1)*Ip(1, 2) + rm(2)*Ip(2, 2) + rm(3)*Ip(3, 2)
197 492418 : cp(3) = rm(1)*Ip(1, 3) + rm(2)*Ip(2, 3) + rm(3)*Ip(3, 3)
198 : ! X Rot
199 492418 : Rot((iparticle - 1)*3 + 1, 1) = (cp(2)*Ip(1, 3) - Ip(1, 2)*cp(3))*mass
200 492418 : Rot((iparticle - 1)*3 + 2, 1) = (cp(2)*Ip(2, 3) - Ip(2, 2)*cp(3))*mass
201 492418 : Rot((iparticle - 1)*3 + 3, 1) = (cp(2)*Ip(3, 3) - Ip(3, 2)*cp(3))*mass
202 : ! Y Rot
203 492418 : Rot((iparticle - 1)*3 + 1, 2) = (cp(3)*Ip(1, 1) - Ip(1, 3)*cp(1))*mass
204 492418 : Rot((iparticle - 1)*3 + 2, 2) = (cp(3)*Ip(2, 1) - Ip(2, 3)*cp(1))*mass
205 492418 : Rot((iparticle - 1)*3 + 3, 2) = (cp(3)*Ip(3, 1) - Ip(3, 3)*cp(1))*mass
206 : ! Z Rot
207 492418 : Rot((iparticle - 1)*3 + 1, 3) = (cp(1)*Ip(1, 2) - Ip(1, 1)*cp(2))*mass
208 492418 : Rot((iparticle - 1)*3 + 2, 3) = (cp(1)*Ip(2, 2) - Ip(2, 1)*cp(2))*mass
209 494583 : 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 8660 : lrot = 1
214 8660 : DO i = 1, 3
215 4438257 : norm = DOT_PRODUCT(Rot(:, i), Rot(:, i))
216 6495 : IF (norm <= thrs_motion) THEN
217 226 : lrot(i) = 0
218 226 : CYCLE
219 : END IF
220 4436711 : Rot(:, i) = Rot(:, i)/SQRT(norm)
221 : ! Clean Rotational modes for spurious/numerical contamination
222 8434 : IF (i < 3) THEN
223 10385 : DO j = 1, i
224 8871269 : 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 7474 : IF (PRESENT(rot_dof)) rot_dof = COUNT(lrot == 1)
230 8668 : dof = dof + COUNT(lrot == 1)
231 2167 : iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
232 2167 : IF (iw > 0) THEN
233 1002 : WRITE (iw, '(T2,A,T71,I10)') 'ROT| Number of rotovibrational vectors', dof
234 1002 : IF (dof == 5) THEN
235 : WRITE (iw, '(T2,A)') &
236 92 : 'ROT| Linear molecule detected'
237 : END IF
238 1002 : IF ((dof == 3) .AND. (.NOT. keep_rotations)) THEN
239 : WRITE (iw, '(T2,A)') &
240 6 : 'ROT| Single atom detected'
241 : END IF
242 : END IF
243 2167 : CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
244 2167 : IF (present_mat) THEN
245 : ! Give back the vectors generating the rototranslating Frame
246 1592 : ALLOCATE (mat(natoms*3, dof))
247 398 : iseq = 0
248 1592 : DO i = 1, 3
249 18366 : mat(:, i) = Tr(:, i)
250 1592 : IF (lrot(i) == 1) THEN
251 1138 : iseq = iseq + 1
252 17956 : mat(:, 3 + iseq) = Rot(:, i)
253 : END IF
254 : END DO
255 : END IF
256 2167 : DEALLOCATE (Tr)
257 2167 : DEALLOCATE (Rot)
258 2167 : CALL timestop(handle)
259 :
260 2167 : 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 213232 : 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 213232 : INTEGER, POINTER :: force_mixing_indices(:), &
304 213232 : force_mixing_labels(:)
305 : LOGICAL :: charge_beta, charge_extended, &
306 : charge_occup, explicit, &
307 : my_extended_xmol_title, new_file, &
308 : print_kind
309 213232 : 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 213232 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
316 : TYPE(section_vals_type), POINTER :: force_env_section, &
317 : force_mixing_restart_section
318 :
319 213232 : CALL timeset(routineN, handle)
320 :
321 213232 : NULLIFY (logger, cell, subsys, my_particles, particle_set)
322 213232 : logger => cp_get_default_logger()
323 213232 : id_label = logger%iter_info%level_name(logger%iter_info%n_rlevel)
324 213232 : my_pos = "APPEND"
325 213232 : my_act = "WRITE"
326 213232 : my_middle = "pos"
327 213232 : my_pk_name = "TRAJECTORY"
328 213232 : IF (PRESENT(middle_name)) my_middle = middle_name
329 213232 : IF (PRESENT(pos)) my_pos = pos
330 213232 : IF (PRESENT(act)) my_act = act
331 213232 : IF (PRESENT(pk_name)) my_pk_name = pk_name
332 :
333 275343 : SELECT CASE (TRIM(my_pk_name))
334 : CASE ("TRAJECTORY", "SHELL_TRAJECTORY", "CORE_TRAJECTORY")
335 62111 : id_dcd = "CORD"
336 62111 : id_wpc = "POS"
337 62111 : id_extxyz = "pos"
338 : CASE ("VELOCITIES", "SHELL_VELOCITIES", "CORE_VELOCITIES")
339 46891 : id_dcd = "VEL "
340 46891 : id_wpc = "VEL"
341 46891 : id_extxyz = "velo"
342 : CASE ("FORCES", "SHELL_FORCES", "CORE_FORCES")
343 62111 : id_dcd = "FRC "
344 62111 : id_wpc = "FORCE"
345 62111 : id_extxyz = "force"
346 : CASE ("FORCE_MIXING_LABELS")
347 42119 : id_dcd = "FML "
348 42119 : id_wpc = "FORCE_MIXING_LABELS"
349 42119 : 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 213232 : "<"//TRIM(my_pk_name)//">")
357 : END SELECT
358 :
359 213232 : charge_occup = .FALSE.
360 213232 : charge_beta = .FALSE.
361 213232 : charge_extended = .FALSE.
362 213232 : print_kind = .FALSE.
363 :
364 213232 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
365 213232 : IF (PRESENT(particles)) THEN
366 27940 : CPASSERT(ASSOCIATED(particles))
367 27940 : my_particles => particles
368 : ELSE
369 185292 : CALL cp_subsys_get(subsys=subsys, particles=my_particles)
370 : END IF
371 213232 : particle_set => my_particles%els
372 213232 : nat = my_particles%n_els
373 213232 : pbc_str = "F F F"
374 213232 : IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
375 213232 : IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
376 213232 : IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
377 :
378 : ! Gather units of measure for output (if available)
379 213232 : 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 171113 : c_val=unit_str)
382 171113 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
383 : END IF
384 :
385 : ! Get the output format
386 213232 : 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 213232 : file_form=my_form, middle_name=TRIM(my_middle), is_new_file=new_file)
390 213232 : IF (traj_unit > 0) THEN
391 : CALL section_vals_val_get(root_section, "MOTION%PRINT%"//TRIM(my_pk_name)//"%FORMAT", &
392 22361 : i_val=outformat)
393 22361 : title = ""
394 4 : SELECT CASE (outformat)
395 : CASE (dump_dcd, dump_dcd_aligned_cell)
396 4 : IF (new_file) THEN
397 : !Lets write the header for the coordinate dcd
398 1 : section_ref = "MOTION%PRINT%"//TRIM(my_pk_name)//"%EACH%"//TRIM(id_label)
399 1 : iskip = section_get_ival(root_section, TRIM(section_ref))
400 : i = INDEX(cp2k_version, "(") - 1
401 : IF (i == -1) i = LEN(cp2k_version)
402 1 : CALL m_timestamp(timestamp)
403 1 : WRITE (UNIT=traj_unit) id_dcd, 0, it, iskip, 0, 0, 0, 0, 0, 0, REAL(dtime, KIND=sp), &
404 2 : 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 1 : " (revision "//TRIM(compile_revision)//")"
407 1 : remark(2) = "REMARK "//TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//timestamp(:19)
408 1 : WRITE (UNIT=traj_unit) SIZE(remark), remark(:)
409 1 : WRITE (UNIT=traj_unit) nat
410 1 : CALL m_flush(traj_unit)
411 : END IF
412 : CASE (dump_xmol)
413 22173 : my_extended_xmol_title = .FALSE.
414 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
415 22173 : l_val=print_kind)
416 22173 : IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
417 : ! This information can be digested by Molden
418 18596 : IF (my_extended_xmol_title) THEN
419 : WRITE (UNIT=title, FMT="(A,I8,A,F12.3,A,F20.10)") &
420 18596 : " i = ", it, ", time = ", time, ", E = ", etot
421 : ELSE
422 3577 : WRITE (UNIT=title, FMT="(A,I8,A,F20.10)") " i = ", it, ", E = ", etot
423 : END IF
424 : CASE (dump_extxyz)
425 115 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", l_val=print_kind)
426 1150 : WRITE (UNIT=cell_str, FMT="(9(1X,F19.10))") cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom
427 115 : WRITE (UNIT=step_str, FMT="(I8)") it
428 115 : WRITE (UNIT=time_str, FMT="(F12.3)") time
429 115 : 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 115 : ' Energy='//TRIM(ADJUSTL(etot_str))
437 : CASE (dump_atomic)
438 : ! Do nothing
439 : CASE (dump_pdb)
440 59 : IF (id_wpc == "POS") THEN
441 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_OCCUP", &
442 59 : l_val=charge_occup)
443 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_BETA", &
444 59 : l_val=charge_beta)
445 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_EXTENDED", &
446 59 : l_val=charge_extended)
447 236 : i = COUNT([charge_occup, charge_beta, charge_extended])
448 59 : IF (i > 1) &
449 0 : CPABORT("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected, ")
450 : END IF
451 59 : IF (new_file) THEN
452 4 : CALL m_timestamp(timestamp)
453 : ! COLUMNS DATA TYPE FIELD DEFINITION
454 : ! 1 - 6 Record name "TITLE "
455 : ! 9 - 10 Continuation continuation Allows concatenation
456 : ! 11 - 70 String title Title of the experiment
457 : WRITE (UNIT=traj_unit, FMT="(A6,T11,A)") &
458 4 : "TITLE ", "PDB file created by "//TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")", &
459 8 : "AUTHOR", TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//timestamp(:19)
460 : END IF
461 59 : my_extended_xmol_title = .FALSE.
462 59 : IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
463 0 : IF (my_extended_xmol_title) THEN
464 : WRITE (UNIT=title, FMT="(A,I0,A,F0.3,A,F0.10)") &
465 0 : "Step ", it, ", time = ", time, ", E = ", etot
466 : ELSE
467 : WRITE (UNIT=title, FMT="(A,I0,A,F0.10)") &
468 59 : "Step ", it, ", E = ", etot
469 : END IF
470 : CASE DEFAULT
471 22361 : CPABORT("Unknown output format")
472 : END SELECT
473 22361 : IF (TRIM(my_pk_name) == "FORCE_MIXING_LABELS") THEN
474 453 : ALLOCATE (fml_array(3*SIZE(particle_set)))
475 151 : fml_array = 0.0_dp
476 151 : CALL force_env_get(force_env, force_env_section=force_env_section)
477 : force_mixing_restart_section => section_vals_get_subs_vals(force_env_section, &
478 : "QMMM%FORCE_MIXING%RESTART_INFO", &
479 151 : can_return_null=.TRUE.)
480 151 : IF (ASSOCIATED(force_mixing_restart_section)) THEN
481 151 : CALL section_vals_get(force_mixing_restart_section, explicit=explicit)
482 151 : IF (explicit) THEN
483 0 : CALL section_vals_val_get(force_mixing_restart_section, "INDICES", i_vals=force_mixing_indices)
484 0 : CALL section_vals_val_get(force_mixing_restart_section, "LABELS", i_vals=force_mixing_labels)
485 0 : DO i = 1, SIZE(force_mixing_indices)
486 0 : ii = force_mixing_indices(i)
487 0 : CPASSERT(ii <= SIZE(particle_set))
488 0 : fml_array((ii - 1)*3 + 1:(ii - 1)*3 + 3) = force_mixing_labels(i)
489 : END DO
490 : END IF
491 : END IF
492 : CALL write_particle_coordinates(particle_set, traj_unit, outformat, TRIM(id_wpc), TRIM(title), cell, &
493 151 : array=fml_array, print_kind=print_kind)
494 151 : DEALLOCATE (fml_array)
495 : ELSE
496 : CALL write_particle_coordinates(particle_set, traj_unit, outformat, TRIM(id_wpc), TRIM(title), cell, &
497 : unit_conv=unit_conv, print_kind=print_kind, &
498 : charge_occup=charge_occup, &
499 : charge_beta=charge_beta, &
500 22210 : charge_extended=charge_extended)
501 : END IF
502 : END IF
503 :
504 213232 : CALL cp_print_key_finished_output(traj_unit, logger, root_section, "MOTION%PRINT%"//TRIM(my_pk_name))
505 :
506 213232 : CALL timestop(handle)
507 :
508 426464 : END SUBROUTINE write_trajectory
509 :
510 : ! **************************************************************************************************
511 : !> \brief Info on the unit to be opened to dump MD informations
512 : !> \param section ...
513 : !> \param path ...
514 : !> \param my_form ...
515 : !> \param my_ext ...
516 : !> \author Teodoro Laino - University of Zurich - 07.2007
517 : ! **************************************************************************************************
518 213270 : SUBROUTINE get_output_format(section, path, my_form, my_ext)
519 :
520 : TYPE(section_vals_type), POINTER :: section
521 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: path
522 : CHARACTER(LEN=*), INTENT(OUT) :: my_form, my_ext
523 :
524 : INTEGER :: output_format
525 :
526 213308 : IF (PRESENT(path)) THEN
527 213232 : CALL section_vals_val_get(section, TRIM(path)//"%FORMAT", i_val=output_format)
528 : ELSE
529 38 : CALL section_vals_val_get(section, "FORMAT", i_val=output_format)
530 : END IF
531 :
532 14 : SELECT CASE (output_format)
533 : CASE (dump_dcd, dump_dcd_aligned_cell)
534 14 : my_form = "UNFORMATTED"
535 14 : my_ext = ".dcd"
536 : CASE (dump_pdb)
537 118 : my_form = "FORMATTED"
538 118 : my_ext = ".pdb"
539 : CASE DEFAULT
540 213138 : my_form = "FORMATTED"
541 426408 : my_ext = ".xyz"
542 : END SELECT
543 :
544 213270 : END SUBROUTINE get_output_format
545 :
546 : ! **************************************************************************************************
547 : !> \brief Prints the Stress Tensor
548 : !> \param virial ...
549 : !> \param cell ...
550 : !> \param motion_section ...
551 : !> \param itimes ...
552 : !> \param time ...
553 : !> \param pos ...
554 : !> \param act ...
555 : !> \date 02.2008
556 : !> \author Teodoro Laino [tlaino] - University of Zurich
557 : !> \version 1.0
558 : ! **************************************************************************************************
559 50169 : SUBROUTINE write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
560 :
561 : TYPE(virial_type), POINTER :: virial
562 : TYPE(cell_type), POINTER :: cell
563 : TYPE(section_vals_type), POINTER :: motion_section
564 : INTEGER, INTENT(IN) :: itimes
565 : REAL(KIND=dp), INTENT(IN) :: time
566 : CHARACTER(LEN=default_string_length), INTENT(IN), &
567 : OPTIONAL :: pos, act
568 :
569 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
570 : INTEGER :: output_unit
571 : LOGICAL :: new_file
572 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_total_bar
573 : TYPE(cp_logger_type), POINTER :: logger
574 :
575 50169 : NULLIFY (logger)
576 50169 : logger => cp_get_default_logger()
577 :
578 50169 : IF (virial%pv_availability) THEN
579 8778 : my_pos = "APPEND"
580 8778 : my_act = "WRITE"
581 8778 : IF (PRESENT(pos)) my_pos = pos
582 8778 : IF (PRESENT(act)) my_act = act
583 : output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%STRESS", &
584 : extension=".stress", file_position=my_pos, &
585 : file_action=my_act, file_form="FORMATTED", &
586 8778 : is_new_file=new_file)
587 : ELSE
588 41391 : output_unit = 0
589 : END IF
590 :
591 50169 : IF (output_unit > 0) THEN
592 1284 : IF (new_file) THEN
593 : WRITE (UNIT=output_unit, FMT='(A,9(12X,A2," [bar]"),6X,A)') &
594 64 : "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
595 : END IF
596 1284 : pv_total_bar(1, 1) = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
597 1284 : pv_total_bar(1, 2) = cp_unit_from_cp2k(virial%pv_total(1, 2)/cell%deth, "bar")
598 1284 : pv_total_bar(1, 3) = cp_unit_from_cp2k(virial%pv_total(1, 3)/cell%deth, "bar")
599 1284 : pv_total_bar(2, 1) = cp_unit_from_cp2k(virial%pv_total(2, 1)/cell%deth, "bar")
600 1284 : pv_total_bar(2, 2) = cp_unit_from_cp2k(virial%pv_total(2, 2)/cell%deth, "bar")
601 1284 : pv_total_bar(2, 3) = cp_unit_from_cp2k(virial%pv_total(2, 3)/cell%deth, "bar")
602 1284 : pv_total_bar(3, 1) = cp_unit_from_cp2k(virial%pv_total(3, 1)/cell%deth, "bar")
603 1284 : pv_total_bar(3, 2) = cp_unit_from_cp2k(virial%pv_total(3, 2)/cell%deth, "bar")
604 1284 : pv_total_bar(3, 3) = cp_unit_from_cp2k(virial%pv_total(3, 3)/cell%deth, "bar")
605 1284 : WRITE (UNIT=output_unit, FMT='(I8,F12.3,9(1X,F19.10))') itimes, time, &
606 1284 : pv_total_bar(1, 1), pv_total_bar(1, 2), pv_total_bar(1, 3), &
607 1284 : pv_total_bar(2, 1), pv_total_bar(2, 2), pv_total_bar(2, 3), &
608 2568 : pv_total_bar(3, 1), pv_total_bar(3, 2), pv_total_bar(3, 3)
609 1284 : CALL m_flush(output_unit)
610 : END IF
611 :
612 50169 : IF (virial%pv_availability) THEN
613 : CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
614 8778 : "PRINT%STRESS")
615 : END IF
616 :
617 50169 : END SUBROUTINE write_stress_tensor_to_file
618 :
619 : ! **************************************************************************************************
620 : !> \brief Prints the Simulation Cell
621 : !> \param cell ...
622 : !> \param motion_section ...
623 : !> \param itimes ...
624 : !> \param time ...
625 : !> \param pos ...
626 : !> \param act ...
627 : !> \date 02.2008
628 : !> \author Teodoro Laino [tlaino] - University of Zurich
629 : !> \version 1.0
630 : ! **************************************************************************************************
631 50169 : SUBROUTINE write_simulation_cell(cell, motion_section, itimes, time, pos, act)
632 :
633 : TYPE(cell_type), POINTER :: cell
634 : TYPE(section_vals_type), POINTER :: motion_section
635 : INTEGER, INTENT(IN) :: itimes
636 : REAL(KIND=dp), INTENT(IN) :: time
637 : CHARACTER(LEN=default_string_length), INTENT(IN), &
638 : OPTIONAL :: pos, act
639 :
640 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
641 : INTEGER :: output_unit
642 : LOGICAL :: new_file
643 : TYPE(cp_logger_type), POINTER :: logger
644 :
645 50169 : NULLIFY (logger)
646 50169 : logger => cp_get_default_logger()
647 :
648 50169 : my_pos = "APPEND"
649 50169 : my_act = "WRITE"
650 50169 : IF (PRESENT(pos)) my_pos = pos
651 50169 : IF (PRESENT(act)) my_act = act
652 :
653 : output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%CELL", &
654 : extension=".cell", file_position=my_pos, &
655 : file_action=my_act, file_form="FORMATTED", &
656 50169 : is_new_file=new_file)
657 :
658 50169 : IF (output_unit > 0) THEN
659 2095 : IF (new_file) THEN
660 : WRITE (UNIT=output_unit, FMT='(A,9(7X,A2," [Angstrom]"),6X,A)') &
661 68 : "# Step Time [fs]", "Ax", "Ay", "Az", "Bx", "By", "Bz", "Cx", "Cy", "Cz", &
662 136 : "Volume [Angstrom^3]"
663 : END IF
664 2095 : WRITE (UNIT=output_unit, FMT="(I8,F12.3,9(1X,F19.10),1X,F24.10)") itimes, time, &
665 2095 : cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom, &
666 2095 : cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom, &
667 2095 : cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom, &
668 4190 : cell%deth*angstrom*angstrom*angstrom
669 2095 : CALL m_flush(output_unit)
670 : END IF
671 :
672 : CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
673 50169 : "PRINT%CELL")
674 :
675 50169 : END SUBROUTINE write_simulation_cell
676 :
677 : END MODULE motion_utils
|
