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