Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 2145 : 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 2145 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rot, Tr
100 : TYPE(cp_logger_type), POINTER :: logger
101 :
102 2145 : CALL timeset(routineN, handle)
103 2145 : logger => cp_get_default_logger()
104 2145 : present_mat = PRESENT(mat)
105 2145 : CPASSERT(ASSOCIATED(particles))
106 2145 : IF (present_mat) THEN
107 376 : CPASSERT(.NOT. ASSOCIATED(mat))
108 : END IF
109 2145 : IF (.NOT. keep_rotations) THEN
110 2143 : rcom = 0.0_dp
111 2143 : masst = 0.0_dp
112 : ! Center of mass
113 489875 : DO iparticle = 1, natoms
114 487732 : mass = 1.0_dp
115 487732 : IF (mass_weighted) mass = particles(iparticle)%atomic_kind%mass
116 487732 : CPASSERT(mass >= 0.0_dp)
117 487732 : masst = masst + mass
118 1953071 : rcom = particles(iparticle)%r*mass + rcom
119 : END DO
120 2143 : CPASSERT(masst > 0.0_dp)
121 8572 : rcom = rcom/masst
122 : ! Intertia Tensor
123 2143 : Ip = 0.0_dp
124 489875 : DO iparticle = 1, natoms
125 487732 : mass = 1.0_dp
126 487732 : IF (mass_weighted) mass = particles(iparticle)%atomic_kind%mass
127 1950928 : rm = particles(iparticle)%r - rcom
128 487732 : Ip(1, 1) = Ip(1, 1) + mass*(rm(2)**2 + rm(3)**2)
129 487732 : Ip(2, 2) = Ip(2, 2) + mass*(rm(1)**2 + rm(3)**2)
130 487732 : Ip(3, 3) = Ip(3, 3) + mass*(rm(1)**2 + rm(2)**2)
131 487732 : Ip(1, 2) = Ip(1, 2) - mass*(rm(1)*rm(2))
132 487732 : Ip(1, 3) = Ip(1, 3) - mass*(rm(1)*rm(3))
133 489875 : Ip(2, 3) = Ip(2, 3) - mass*(rm(2)*rm(3))
134 : END DO
135 : ! Diagonalize the Inertia Tensor
136 2143 : CALL diamat_all(Ip, Ip_eigval)
137 2143 : IF (PRESENT(inertia)) inertia = Ip_eigval
138 2143 : iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
139 2143 : IF (iw > 0) THEN
140 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
141 992 : 'ROT| Rotational analysis information'
142 : WRITE (UNIT=iw, FMT='(T2,A)') &
143 992 : 'ROT| Principal axes and moments of inertia [a.u.]'
144 : WRITE (UNIT=iw, FMT='(T2,A,T14,3(1X,I19))') &
145 992 : 'ROT|', 1, 2, 3
146 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
147 992 : 'ROT| Eigenvalues', Ip_eigval(1:3)
148 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
149 992 : 'ROT| x', Ip(1, 1:3)
150 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
151 992 : 'ROT| y', Ip(2, 1:3)
152 : WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
153 992 : 'ROT| z', Ip(3, 1:3)
154 : END IF
155 2143 : CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
156 2143 : iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO/COORDINATES", extension=".vibLog")
157 2143 : 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 2143 : CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO/COORDINATES")
166 : END IF
167 : ! Build up the Translational vectors
168 8580 : ALLOCATE (Tr(natoms*3, 3))
169 4398222 : Tr = 0.0_dp
170 8580 : DO k = 1, 3
171 : iseq = 0
172 1471794 : DO iparticle = 1, natoms
173 1463214 : mass = 1.0_dp
174 1463214 : IF (mass_weighted) mass = SQRT(particles(iparticle)%atomic_kind%mass)
175 5859291 : DO j = 1, 3
176 4389642 : iseq = iseq + 1
177 5852856 : IF (j == k) Tr(iseq, k) = mass
178 : END DO
179 : END DO
180 : END DO
181 : ! Normalize Translations
182 8580 : DO i = 1, 3
183 4396077 : norm = SQRT(DOT_PRODUCT(Tr(:, i), Tr(:, i)))
184 4398222 : Tr(:, i) = Tr(:, i)/norm
185 : END DO
186 2145 : dof = 3
187 : ! Build up the Rotational vectors
188 4290 : ALLOCATE (Rot(natoms*3, 3))
189 2145 : lrot = 0
190 2145 : IF (.NOT. keep_rotations) THEN
191 489875 : DO iparticle = 1, natoms
192 487732 : mass = 1.0_dp
193 487732 : IF (mass_weighted) mass = SQRT(particles(iparticle)%atomic_kind%mass)
194 1950928 : rm = particles(iparticle)%r - rcom
195 487732 : cp(1) = rm(1)*Ip(1, 1) + rm(2)*Ip(2, 1) + rm(3)*Ip(3, 1)
196 487732 : cp(2) = rm(1)*Ip(1, 2) + rm(2)*Ip(2, 2) + rm(3)*Ip(3, 2)
197 487732 : cp(3) = rm(1)*Ip(1, 3) + rm(2)*Ip(2, 3) + rm(3)*Ip(3, 3)
198 : ! X Rot
199 487732 : Rot((iparticle - 1)*3 + 1, 1) = (cp(2)*Ip(1, 3) - Ip(1, 2)*cp(3))*mass
200 487732 : Rot((iparticle - 1)*3 + 2, 1) = (cp(2)*Ip(2, 3) - Ip(2, 2)*cp(3))*mass
201 487732 : Rot((iparticle - 1)*3 + 3, 1) = (cp(2)*Ip(3, 3) - Ip(3, 2)*cp(3))*mass
202 : ! Y Rot
203 487732 : Rot((iparticle - 1)*3 + 1, 2) = (cp(3)*Ip(1, 1) - Ip(1, 3)*cp(1))*mass
204 487732 : Rot((iparticle - 1)*3 + 2, 2) = (cp(3)*Ip(2, 1) - Ip(2, 3)*cp(1))*mass
205 487732 : Rot((iparticle - 1)*3 + 3, 2) = (cp(3)*Ip(3, 1) - Ip(3, 3)*cp(1))*mass
206 : ! Z Rot
207 487732 : Rot((iparticle - 1)*3 + 1, 3) = (cp(1)*Ip(1, 2) - Ip(1, 1)*cp(2))*mass
208 487732 : Rot((iparticle - 1)*3 + 2, 3) = (cp(1)*Ip(2, 2) - Ip(2, 1)*cp(2))*mass
209 489875 : 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 8572 : lrot = 1
214 8572 : DO i = 1, 3
215 4396017 : norm = DOT_PRODUCT(Rot(:, i), Rot(:, i))
216 6429 : IF (norm <= thrs_motion) THEN
217 226 : lrot(i) = 0
218 226 : CYCLE
219 : END IF
220 4394471 : Rot(:, i) = Rot(:, i)/SQRT(norm)
221 : ! Clean Rotational modes for spurious/numerical contamination
222 8346 : IF (i < 3) THEN
223 10275 : DO j = 1, i
224 8786811 : 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 7452 : IF (PRESENT(rot_dof)) rot_dof = COUNT(lrot == 1)
230 8580 : dof = dof + COUNT(lrot == 1)
231 2145 : iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
232 2145 : IF (iw > 0) THEN
233 992 : WRITE (iw, '(T2,A,T71,I10)') 'ROT| Number of rotovibrational vectors', dof
234 992 : IF (dof == 5) THEN
235 : WRITE (iw, '(T2,A)') &
236 92 : 'ROT| Linear molecule detected'
237 : END IF
238 992 : 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 2145 : CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
244 2145 : IF (present_mat) THEN
245 : ! Give back the vectors generating the rototranslating Frame
246 1504 : ALLOCATE (mat(natoms*3, dof))
247 376 : iseq = 0
248 1504 : DO i = 1, 3
249 17310 : mat(:, i) = Tr(:, i)
250 1504 : IF (lrot(i) == 1) THEN
251 1072 : iseq = iseq + 1
252 16900 : mat(:, 3 + iseq) = Rot(:, i)
253 : END IF
254 : END DO
255 : END IF
256 2145 : DEALLOCATE (Tr)
257 2145 : DEALLOCATE (Rot)
258 2145 : CALL timestop(handle)
259 :
260 2145 : 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 224850 : SUBROUTINE write_trajectory(force_env, root_section, it, time, dtime, etot, pk_name, &
281 : pos, act, middle_name, particles, extended_xmol_title)
282 : TYPE(force_env_type), POINTER :: force_env
283 : TYPE(section_vals_type), POINTER :: root_section
284 : INTEGER, INTENT(IN) :: it
285 : REAL(KIND=dp), INTENT(IN) :: time, dtime, etot
286 : CHARACTER(LEN=*), OPTIONAL :: pk_name
287 : CHARACTER(LEN=default_string_length), OPTIONAL :: pos, act
288 : CHARACTER(LEN=*), OPTIONAL :: middle_name
289 : TYPE(particle_list_type), OPTIONAL, POINTER :: particles
290 : LOGICAL, INTENT(IN), OPTIONAL :: extended_xmol_title
291 :
292 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_trajectory'
293 :
294 : CHARACTER(LEN=1024) :: cell_str, title
295 : CHARACTER(LEN=4) :: id_dcd
296 : CHARACTER(LEN=80), DIMENSION(2) :: remark
297 : CHARACTER(LEN=default_string_length) :: etot_str, id_extxyz, id_label, id_wpc, my_act, &
298 : my_ext, my_form, my_middle, my_pk_name, my_pos, section_ref, step_str, time_str, unit_str
299 : CHARACTER(LEN=timestamp_length) :: timestamp
300 : INTEGER :: handle, i, ii, iskip, nat, outformat, &
301 : traj_unit
302 224850 : INTEGER, POINTER :: force_mixing_indices(:), &
303 224850 : force_mixing_labels(:)
304 : LOGICAL :: charge_beta, charge_extended, &
305 : charge_occup, explicit, &
306 : my_extended_xmol_title, new_file, &
307 : print_kind
308 224850 : REAL(dp), ALLOCATABLE :: fml_array(:)
309 : REAL(KIND=dp) :: unit_conv
310 : TYPE(cell_type), POINTER :: cell
311 : TYPE(cp_logger_type), POINTER :: logger
312 : TYPE(cp_subsys_type), POINTER :: subsys
313 : TYPE(particle_list_type), POINTER :: my_particles
314 224850 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
315 : TYPE(section_vals_type), POINTER :: force_env_section, &
316 : force_mixing_restart_section
317 :
318 224850 : CALL timeset(routineN, handle)
319 :
320 224850 : NULLIFY (logger, cell, subsys, my_particles, particle_set)
321 224850 : logger => cp_get_default_logger()
322 224850 : id_label = logger%iter_info%level_name(logger%iter_info%n_rlevel)
323 224850 : my_pos = "APPEND"
324 224850 : my_act = "WRITE"
325 224850 : my_middle = "pos"
326 224850 : my_pk_name = "TRAJECTORY"
327 224850 : IF (PRESENT(middle_name)) my_middle = middle_name
328 224850 : IF (PRESENT(pos)) my_pos = pos
329 224850 : IF (PRESENT(act)) my_act = act
330 224850 : IF (PRESENT(pk_name)) my_pk_name = pk_name
331 :
332 292740 : SELECT CASE (TRIM(my_pk_name))
333 : CASE ("TRAJECTORY", "SHELL_TRAJECTORY", "CORE_TRAJECTORY")
334 67890 : id_dcd = "CORD"
335 67890 : id_wpc = "POS"
336 67890 : id_extxyz = "pos"
337 : CASE ("VELOCITIES", "SHELL_VELOCITIES", "CORE_VELOCITIES")
338 46921 : id_dcd = "VEL "
339 46921 : id_wpc = "VEL"
340 46921 : id_extxyz = "velo"
341 : CASE ("FORCES", "SHELL_FORCES", "CORE_FORCES")
342 67890 : id_dcd = "FRC "
343 67890 : id_wpc = "FORCE"
344 67890 : id_extxyz = "force"
345 : CASE ("FORCE_MIXING_LABELS")
346 42149 : id_dcd = "FML "
347 42149 : id_wpc = "FORCE_MIXING_LABELS"
348 42149 : id_extxyz = "force_mixing_label" ! non-standard
349 : CASE DEFAULT
350 224850 : CPABORT("")
351 : END SELECT
352 :
353 224850 : charge_occup = .FALSE.
354 224850 : charge_beta = .FALSE.
355 224850 : charge_extended = .FALSE.
356 224850 : print_kind = .FALSE.
357 :
358 224850 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
359 224850 : IF (PRESENT(particles)) THEN
360 30828 : CPASSERT(ASSOCIATED(particles))
361 30828 : my_particles => particles
362 : ELSE
363 194022 : CALL cp_subsys_get(subsys=subsys, particles=my_particles)
364 : END IF
365 224850 : particle_set => my_particles%els
366 224850 : nat = my_particles%n_els
367 :
368 : ! Gather units of measure for output (if available)
369 224850 : IF (TRIM(my_pk_name) /= "FORCE_MIXING_LABELS") THEN
370 : CALL section_vals_val_get(root_section, "MOTION%PRINT%"//TRIM(my_pk_name)//"%UNIT", &
371 182701 : c_val=unit_str)
372 182701 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
373 : END IF
374 :
375 : ! Get the output format
376 224850 : CALL get_output_format(root_section, "MOTION%PRINT%"//TRIM(my_pk_name), my_form, my_ext)
377 : traj_unit = cp_print_key_unit_nr(logger, root_section, "MOTION%PRINT%"//TRIM(my_pk_name), &
378 : extension=my_ext, file_position=my_pos, file_action=my_act, &
379 224850 : file_form=my_form, middle_name=TRIM(my_middle), is_new_file=new_file)
380 224850 : IF (traj_unit > 0) THEN
381 : CALL section_vals_val_get(root_section, "MOTION%PRINT%"//TRIM(my_pk_name)//"%FORMAT", &
382 23037 : i_val=outformat)
383 23037 : title = ""
384 4 : SELECT CASE (outformat)
385 : CASE (dump_dcd, dump_dcd_aligned_cell)
386 4 : IF (new_file) THEN
387 : !Lets write the header for the coordinate dcd
388 1 : section_ref = "MOTION%PRINT%"//TRIM(my_pk_name)//"%EACH%"//TRIM(id_label)
389 1 : iskip = section_get_ival(root_section, TRIM(section_ref))
390 : i = INDEX(cp2k_version, "(") - 1
391 : IF (i == -1) i = LEN(cp2k_version)
392 1 : CALL m_timestamp(timestamp)
393 1 : WRITE (UNIT=traj_unit) id_dcd, 0, it, iskip, 0, 0, 0, 0, 0, 0, REAL(dtime, KIND=sp), &
394 2 : 1, 0, 0, 0, 0, 0, 0, 0, 0, 24
395 : remark(1) = "REMARK "//id_dcd//" DCD file created by "//TRIM(cp2k_version(:i))// &
396 1 : " (revision "//TRIM(compile_revision)//")"
397 1 : remark(2) = "REMARK "//TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//timestamp(:19)
398 1 : WRITE (UNIT=traj_unit) SIZE(remark), remark(:)
399 1 : WRITE (UNIT=traj_unit) nat
400 1 : CALL m_flush(traj_unit)
401 : END IF
402 : CASE (dump_xmol)
403 22960 : my_extended_xmol_title = .FALSE.
404 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
405 22960 : l_val=print_kind)
406 22960 : IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
407 : ! This information can be digested by Molden
408 18619 : IF (my_extended_xmol_title) THEN
409 : WRITE (UNIT=title, FMT="(A,I8,A,F12.3,A,F20.10)") &
410 18619 : " i = ", it, ", time = ", time, ", E = ", etot
411 : ELSE
412 4341 : WRITE (UNIT=title, FMT="(A,I8,A,F20.10)") " i = ", it, ", E = ", etot
413 : END IF
414 : CASE (dump_extxyz)
415 4 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", l_val=print_kind)
416 40 : WRITE (UNIT=cell_str, FMT="(9(1X,F10.3))") cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom
417 4 : WRITE (UNIT=step_str, FMT="(I8)") it
418 4 : WRITE (UNIT=time_str, FMT="(F12.3)") time
419 4 : WRITE (UNIT=etot_str, FMT="(F20.10)") etot
420 : WRITE (UNIT=title, FMT="(A)") &
421 : 'Lattice="'//TRIM(ADJUSTL(cell_str))//'"'// &
422 : ' Properties="species:S:1:'//TRIM(id_extxyz)//':R:3"'// &
423 : ' Step='//TRIM(ADJUSTL(step_str))// &
424 : ' Time='//TRIM(ADJUSTL(time_str))// &
425 4 : ' Energy='//TRIM(ADJUSTL(etot_str))
426 : CASE (dump_atomic)
427 : ! Do nothing
428 : CASE (dump_pdb)
429 59 : IF (id_wpc == "POS") THEN
430 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_OCCUP", &
431 59 : l_val=charge_occup)
432 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_BETA", &
433 59 : l_val=charge_beta)
434 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%CHARGE_EXTENDED", &
435 59 : l_val=charge_extended)
436 236 : i = COUNT([charge_occup, charge_beta, charge_extended])
437 59 : IF (i > 1) &
438 0 : CPABORT("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected, ")
439 : END IF
440 59 : IF (new_file) THEN
441 4 : CALL m_timestamp(timestamp)
442 : ! COLUMNS DATA TYPE FIELD DEFINITION
443 : ! 1 - 6 Record name "TITLE "
444 : ! 9 - 10 Continuation continuation Allows concatenation
445 : ! 11 - 70 String title Title of the experiment
446 : WRITE (UNIT=traj_unit, FMT="(A6,T11,A)") &
447 4 : "TITLE ", "PDB file created by "//TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")", &
448 8 : "AUTHOR", TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//timestamp(:19)
449 : END IF
450 59 : my_extended_xmol_title = .FALSE.
451 59 : IF (PRESENT(extended_xmol_title)) my_extended_xmol_title = extended_xmol_title
452 0 : IF (my_extended_xmol_title) THEN
453 : WRITE (UNIT=title, FMT="(A,I0,A,F0.3,A,F0.10)") &
454 0 : "Step ", it, ", time = ", time, ", E = ", etot
455 : ELSE
456 : WRITE (UNIT=title, FMT="(A,I0,A,F0.10)") &
457 59 : "Step ", it, ", E = ", etot
458 : END IF
459 : CASE DEFAULT
460 23037 : CPABORT("")
461 : END SELECT
462 23037 : IF (TRIM(my_pk_name) == "FORCE_MIXING_LABELS") THEN
463 453 : ALLOCATE (fml_array(3*SIZE(particle_set)))
464 25978 : fml_array = 0.0_dp
465 151 : CALL force_env_get(force_env, force_env_section=force_env_section)
466 : force_mixing_restart_section => section_vals_get_subs_vals(force_env_section, &
467 : "QMMM%FORCE_MIXING%RESTART_INFO", &
468 151 : can_return_null=.TRUE.)
469 151 : IF (ASSOCIATED(force_mixing_restart_section)) THEN
470 151 : CALL section_vals_get(force_mixing_restart_section, explicit=explicit)
471 151 : IF (explicit) THEN
472 0 : CALL section_vals_val_get(force_mixing_restart_section, "INDICES", i_vals=force_mixing_indices)
473 0 : CALL section_vals_val_get(force_mixing_restart_section, "LABELS", i_vals=force_mixing_labels)
474 0 : DO i = 1, SIZE(force_mixing_indices)
475 0 : ii = force_mixing_indices(i)
476 0 : CPASSERT(ii <= SIZE(particle_set))
477 0 : fml_array((ii - 1)*3 + 1:(ii - 1)*3 + 3) = force_mixing_labels(i)
478 : END DO
479 : END IF
480 : END IF
481 : CALL write_particle_coordinates(particle_set, traj_unit, outformat, TRIM(id_wpc), TRIM(title), cell, &
482 151 : array=fml_array, print_kind=print_kind)
483 151 : DEALLOCATE (fml_array)
484 : ELSE
485 : CALL write_particle_coordinates(particle_set, traj_unit, outformat, TRIM(id_wpc), TRIM(title), cell, &
486 : unit_conv=unit_conv, print_kind=print_kind, &
487 : charge_occup=charge_occup, &
488 : charge_beta=charge_beta, &
489 22886 : charge_extended=charge_extended)
490 : END IF
491 : END IF
492 :
493 224850 : CALL cp_print_key_finished_output(traj_unit, logger, root_section, "MOTION%PRINT%"//TRIM(my_pk_name))
494 :
495 224850 : CALL timestop(handle)
496 :
497 449700 : END SUBROUTINE write_trajectory
498 :
499 : ! **************************************************************************************************
500 : !> \brief Info on the unit to be opened to dump MD informations
501 : !> \param section ...
502 : !> \param path ...
503 : !> \param my_form ...
504 : !> \param my_ext ...
505 : !> \author Teodoro Laino - University of Zurich - 07.2007
506 : ! **************************************************************************************************
507 224876 : SUBROUTINE get_output_format(section, path, my_form, my_ext)
508 :
509 : TYPE(section_vals_type), POINTER :: section
510 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: path
511 : CHARACTER(LEN=*), INTENT(OUT) :: my_form, my_ext
512 :
513 : INTEGER :: output_format
514 :
515 224902 : IF (PRESENT(path)) THEN
516 224850 : CALL section_vals_val_get(section, TRIM(path)//"%FORMAT", i_val=output_format)
517 : ELSE
518 26 : CALL section_vals_val_get(section, "FORMAT", i_val=output_format)
519 : END IF
520 :
521 14 : SELECT CASE (output_format)
522 : CASE (dump_dcd, dump_dcd_aligned_cell)
523 14 : my_form = "UNFORMATTED"
524 14 : my_ext = ".dcd"
525 : CASE (dump_pdb)
526 118 : my_form = "FORMATTED"
527 118 : my_ext = ".pdb"
528 : CASE DEFAULT
529 224744 : my_form = "FORMATTED"
530 449620 : my_ext = ".xyz"
531 : END SELECT
532 :
533 224876 : END SUBROUTINE get_output_format
534 :
535 : ! **************************************************************************************************
536 : !> \brief Prints the Stress Tensor
537 : !> \param virial ...
538 : !> \param cell ...
539 : !> \param motion_section ...
540 : !> \param itimes ...
541 : !> \param time ...
542 : !> \param pos ...
543 : !> \param act ...
544 : !> \date 02.2008
545 : !> \author Teodoro Laino [tlaino] - University of Zurich
546 : !> \version 1.0
547 : ! **************************************************************************************************
548 54490 : SUBROUTINE write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
549 :
550 : TYPE(virial_type), POINTER :: virial
551 : TYPE(cell_type), POINTER :: cell
552 : TYPE(section_vals_type), POINTER :: motion_section
553 : INTEGER, INTENT(IN) :: itimes
554 : REAL(KIND=dp), INTENT(IN) :: time
555 : CHARACTER(LEN=default_string_length), INTENT(IN), &
556 : OPTIONAL :: pos, act
557 :
558 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
559 : INTEGER :: output_unit
560 : LOGICAL :: new_file
561 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_total_bar
562 : TYPE(cp_logger_type), POINTER :: logger
563 :
564 54490 : NULLIFY (logger)
565 54490 : logger => cp_get_default_logger()
566 :
567 54490 : IF (virial%pv_availability) THEN
568 13286 : my_pos = "APPEND"
569 13286 : my_act = "WRITE"
570 13286 : IF (PRESENT(pos)) my_pos = pos
571 13286 : IF (PRESENT(act)) my_act = act
572 : output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%STRESS", &
573 : extension=".stress", file_position=my_pos, &
574 : file_action=my_act, file_form="FORMATTED", &
575 13286 : is_new_file=new_file)
576 : ELSE
577 41204 : output_unit = 0
578 : END IF
579 :
580 54490 : IF (output_unit > 0) THEN
581 1490 : IF (new_file) THEN
582 : WRITE (UNIT=output_unit, FMT='(A,9(12X,A2," [bar]"),6X,A)') &
583 74 : "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
584 : END IF
585 1490 : pv_total_bar(1, 1) = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
586 1490 : pv_total_bar(1, 2) = cp_unit_from_cp2k(virial%pv_total(1, 2)/cell%deth, "bar")
587 1490 : pv_total_bar(1, 3) = cp_unit_from_cp2k(virial%pv_total(1, 3)/cell%deth, "bar")
588 1490 : pv_total_bar(2, 1) = cp_unit_from_cp2k(virial%pv_total(2, 1)/cell%deth, "bar")
589 1490 : pv_total_bar(2, 2) = cp_unit_from_cp2k(virial%pv_total(2, 2)/cell%deth, "bar")
590 1490 : pv_total_bar(2, 3) = cp_unit_from_cp2k(virial%pv_total(2, 3)/cell%deth, "bar")
591 1490 : pv_total_bar(3, 1) = cp_unit_from_cp2k(virial%pv_total(3, 1)/cell%deth, "bar")
592 1490 : pv_total_bar(3, 2) = cp_unit_from_cp2k(virial%pv_total(3, 2)/cell%deth, "bar")
593 1490 : pv_total_bar(3, 3) = cp_unit_from_cp2k(virial%pv_total(3, 3)/cell%deth, "bar")
594 1490 : WRITE (UNIT=output_unit, FMT='(I8,F12.3,9(1X,F19.10))') itimes, time, &
595 1490 : pv_total_bar(1, 1), pv_total_bar(1, 2), pv_total_bar(1, 3), &
596 1490 : pv_total_bar(2, 1), pv_total_bar(2, 2), pv_total_bar(2, 3), &
597 2980 : pv_total_bar(3, 1), pv_total_bar(3, 2), pv_total_bar(3, 3)
598 1490 : CALL m_flush(output_unit)
599 : END IF
600 :
601 54490 : IF (virial%pv_availability) THEN
602 : CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
603 13286 : "PRINT%STRESS")
604 : END IF
605 :
606 54490 : END SUBROUTINE write_stress_tensor_to_file
607 :
608 : ! **************************************************************************************************
609 : !> \brief Prints the Simulation Cell
610 : !> \param cell ...
611 : !> \param motion_section ...
612 : !> \param itimes ...
613 : !> \param time ...
614 : !> \param pos ...
615 : !> \param act ...
616 : !> \date 02.2008
617 : !> \author Teodoro Laino [tlaino] - University of Zurich
618 : !> \version 1.0
619 : ! **************************************************************************************************
620 54490 : SUBROUTINE write_simulation_cell(cell, motion_section, itimes, time, pos, act)
621 :
622 : TYPE(cell_type), POINTER :: cell
623 : TYPE(section_vals_type), POINTER :: motion_section
624 : INTEGER, INTENT(IN) :: itimes
625 : REAL(KIND=dp), INTENT(IN) :: time
626 : CHARACTER(LEN=default_string_length), INTENT(IN), &
627 : OPTIONAL :: pos, act
628 :
629 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
630 : INTEGER :: output_unit
631 : LOGICAL :: new_file
632 : TYPE(cp_logger_type), POINTER :: logger
633 :
634 54490 : NULLIFY (logger)
635 54490 : logger => cp_get_default_logger()
636 :
637 54490 : my_pos = "APPEND"
638 54490 : my_act = "WRITE"
639 54490 : IF (PRESENT(pos)) my_pos = pos
640 54490 : IF (PRESENT(act)) my_act = act
641 :
642 : output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%CELL", &
643 : extension=".cell", file_position=my_pos, &
644 : file_action=my_act, file_form="FORMATTED", &
645 54490 : is_new_file=new_file)
646 :
647 54490 : IF (output_unit > 0) THEN
648 2234 : IF (new_file) THEN
649 : WRITE (UNIT=output_unit, FMT='(A,9(7X,A2," [Angstrom]"),6X,A)') &
650 106 : "# Step Time [fs]", "Ax", "Ay", "Az", "Bx", "By", "Bz", "Cx", "Cy", "Cz", &
651 212 : "Volume [Angstrom^3]"
652 : END IF
653 2234 : WRITE (UNIT=output_unit, FMT="(I8,F12.3,9(1X,F19.10),1X,F24.10)") itimes, time, &
654 2234 : cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom, &
655 2234 : cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom, &
656 2234 : cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom, &
657 4468 : cell%deth*angstrom*angstrom*angstrom
658 2234 : CALL m_flush(output_unit)
659 : END IF
660 :
661 : CALL cp_print_key_finished_output(output_unit, logger, motion_section, &
662 54490 : "PRINT%CELL")
663 :
664 54490 : END SUBROUTINE write_simulation_cell
665 :
666 : END MODULE motion_utils
|
