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