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