Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief I/O subroutines for pint_env
10 : !> \author Lukasz Walewski
11 : !> \date 2009-06-04
12 : ! **************************************************************************************************
13 : MODULE pint_io
14 :
15 : USE cell_types, ONLY: cell_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_get_default_unit_nr,&
18 : cp_logger_type
19 : USE cp_output_handling, ONLY: cp_p_file,&
20 : cp_print_key_finished_output,&
21 : cp_print_key_should_output,&
22 : cp_print_key_unit_nr
23 : USE cp_subsys_types, ONLY: cp_subsys_get,&
24 : cp_subsys_type
25 : USE cp_units, ONLY: cp_unit_from_cp2k
26 : USE f77_interface, ONLY: f_env_add_defaults,&
27 : f_env_rm_defaults,&
28 : f_env_type
29 : USE force_env_types, ONLY: force_env_get
30 : USE input_constants, ONLY: dump_atomic,&
31 : dump_dcd,&
32 : dump_dcd_aligned_cell,&
33 : dump_xmol
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type,&
36 : section_vals_val_get
37 : USE kinds, ONLY: default_string_length,&
38 : dp
39 : USE machine, ONLY: m_flush
40 : USE particle_list_types, ONLY: particle_list_type
41 : USE particle_methods, ONLY: write_particle_coordinates
42 : USE pint_public, ONLY: pint_com_pos
43 : USE pint_transformations, ONLY: pint_u2x
44 : USE pint_types, ONLY: e_conserved_id,&
45 : e_kin_thermo_id,&
46 : e_kin_virial_id,&
47 : e_potential_id,&
48 : pint_env_type
49 : #include "../base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
56 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_io'
57 :
58 : PUBLIC :: pint_write_line
59 : PUBLIC :: pint_write_centroids
60 : PUBLIC :: pint_write_trajectory
61 : PUBLIC :: pint_write_com
62 : PUBLIC :: pint_write_ener
63 : PUBLIC :: pint_write_action
64 : PUBLIC :: pint_write_step_info
65 : PUBLIC :: pint_write_rgyr
66 :
67 : CONTAINS
68 :
69 : ! ***************************************************************************
70 : !> \brief Writes out a line of text to the default output unit.
71 : !> \param line ...
72 : !> \date 2009-07-10
73 : !> \author Lukasz Walewski
74 : ! **************************************************************************************************
75 182 : SUBROUTINE pint_write_line(line)
76 :
77 : CHARACTER(len=*), INTENT(IN) :: line
78 :
79 : CHARACTER(len=default_string_length) :: my_label
80 : INTEGER :: unit_nr
81 : TYPE(cp_logger_type), POINTER :: logger
82 :
83 182 : NULLIFY (logger)
84 182 : logger => cp_get_default_logger()
85 182 : my_label = "PINT|"
86 :
87 182 : IF (logger%para_env%is_source()) THEN
88 101 : unit_nr = cp_logger_get_default_unit_nr(logger)
89 101 : WRITE (unit_nr, '(T2,A)') TRIM(my_label)//" "//TRIM(line)
90 : END IF
91 :
92 182 : END SUBROUTINE pint_write_line
93 :
94 : ! ***************************************************************************
95 : !> \brief Write out the trajectory of the centroid (positions and velocities)
96 : !> \param pint_env ...
97 : !> \par History
98 : !> various bug fixes - hforbert
99 : !> 2010-11-25 rewritten, added support for velocity printing,
100 : !> calc of the stddev of the beads turned off [lwalewski]
101 : !> \author fawzi
102 : ! **************************************************************************************************
103 482 : SUBROUTINE pint_write_centroids(pint_env)
104 : TYPE(pint_env_type), INTENT(IN) :: pint_env
105 :
106 : CHARACTER(len=*), PARAMETER :: routineN = 'pint_write_centroids'
107 : INTEGER, PARAMETER :: n_ids = 2, pos_id = 1, vel_id = 2
108 :
109 : CHARACTER(len=default_string_length) :: ext, form, my_middle_name, unit_str
110 : CHARACTER(len=default_string_length), DIMENSION(2) :: content_id, middle_name, sect_path, title
111 : INTEGER :: handle, handle1, iat, ib, id, idim, &
112 : idir, ierr, outformat, should_output, &
113 : unit_nr
114 : LOGICAL :: new_file, print_kind
115 : REAL(kind=dp) :: nb, ss, unit_conv, vv
116 : TYPE(cell_type), POINTER :: cell
117 : TYPE(cp_logger_type), POINTER :: logger
118 : TYPE(cp_subsys_type), POINTER :: subsys
119 : TYPE(f_env_type), POINTER :: f_env
120 : TYPE(particle_list_type), POINTER :: particles
121 : TYPE(section_vals_type), POINTER :: print_key
122 :
123 482 : CALL timeset(routineN, handle1)
124 :
125 482 : sect_path(pos_id) = "MOTION%PINT%PRINT%CENTROID_POS"
126 482 : sect_path(vel_id) = "MOTION%PINT%PRINT%CENTROID_VEL"
127 482 : middle_name(pos_id) = "centroid-pos"
128 482 : middle_name(vel_id) = "centroid-vel"
129 482 : content_id(pos_id) = "POS"
130 482 : content_id(vel_id) = "VEL"
131 : WRITE (UNIT=title(pos_id), FMT="(A,I8,A,F20.10)") &
132 482 : " i =", pint_env%iter, &
133 4260 : ", E =", SUM(pint_env%e_pot_bead)*pint_env%propagator%physpotscale
134 : WRITE (UNIT=title(vel_id), FMT="(A,I8,A,F20.10,A,F20.10)") &
135 482 : " i =", pint_env%iter, &
136 482 : ", E_trm =", pint_env%energy(e_kin_thermo_id), &
137 964 : ", E_vir =", pint_env%energy(e_kin_virial_id)
138 :
139 482 : NULLIFY (logger)
140 482 : logger => cp_get_default_logger()
141 :
142 482 : CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
143 :
144 : ! iterate over the properties that we know how to print
145 : ! (currently positions and velocities)
146 1446 : DO id = 1, n_ids
147 :
148 : print_key => section_vals_get_subs_vals(pint_env%input, &
149 964 : TRIM(sect_path(id)))
150 :
151 : should_output = cp_print_key_should_output( &
152 : iteration_info=logger%iter_info, &
153 964 : basis_section=print_key)
154 : IF (.NOT. BTEST(should_output, cp_p_file)) CONTINUE
155 :
156 964 : print_kind = .FALSE.
157 :
158 : ! get units of measure for output (if available)
159 : CALL section_vals_val_get(print_key, "UNIT", &
160 964 : c_val=unit_str)
161 964 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
162 :
163 : ! get the format for output
164 964 : CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
165 :
166 0 : SELECT CASE (outformat)
167 : CASE (dump_dcd, dump_dcd_aligned_cell)
168 0 : form = "UNFORMATTED"
169 0 : ext = ".dcd"
170 : CASE (dump_atomic)
171 0 : form = "FORMATTED"
172 0 : ext = ""
173 : CASE (dump_xmol)
174 : CALL section_vals_val_get(print_key, "PRINT_ATOM_KIND", &
175 964 : l_val=print_kind)
176 964 : form = "FORMATTED"
177 964 : ext = ".xyz"
178 : CASE default
179 964 : CPABORT("")
180 : END SELECT
181 :
182 964 : NULLIFY (f_env, cell, subsys)
183 : CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
184 964 : f_env=f_env, handle=handle)
185 : CALL force_env_get(force_env=f_env%force_env, &
186 964 : cell=cell, subsys=subsys)
187 964 : CALL cp_subsys_get(subsys, particles=particles)
188 :
189 : ! calculate and copy the requested property
190 : ! to the particles structure
191 964 : nb = REAL(pint_env%p, dp)
192 964 : idim = 0
193 224740 : DO iat = 1, pint_env%ndim/3
194 896068 : DO idir = 1, 3
195 671328 : idim = idim + 1
196 671328 : ss = 0.0_dp
197 671328 : vv = 0.0_dp
198 : ! ss2=0.0_dp
199 4045104 : DO ib = 1, pint_env%p
200 3373776 : ss = ss + pint_env%x(ib, idim)
201 4045104 : vv = vv + pint_env%v(ib, idim)
202 : ! ss2=ss2+pint_env%x(ib,idim)**2
203 : END DO
204 671328 : particles%els(iat)%r(idir) = ss/nb
205 895104 : particles%els(iat)%v(idir) = vv/nb
206 : ! particles%els(iat)%v(idir)=SQRT(ss2/nb-(ss/nb)**2)
207 : END DO
208 : END DO
209 :
210 : ! set up the output unit number and file name
211 : ! for the current property
212 964 : my_middle_name = TRIM(middle_name(id))
213 : unit_nr = cp_print_key_unit_nr(logger=logger, &
214 : basis_section=print_key, print_key_path="", &
215 : extension=TRIM(ext), middle_name=TRIM(my_middle_name), &
216 964 : local=.FALSE., file_form=form, is_new_file=new_file)
217 :
218 : ! don't write the 0-th frame if the file already exists
219 964 : IF (.NOT. new_file .AND. (pint_env%iter .LE. pint_env%first_step)) THEN
220 : CALL cp_print_key_finished_output(unit_nr, logger, &
221 68 : print_key)
222 : CONTINUE
223 : END IF
224 :
225 : ! actually perform the i/o - on the ionode only
226 964 : IF (unit_nr > 0) THEN
227 :
228 : CALL write_particle_coordinates( &
229 : particles%els, &
230 : iunit=unit_nr, &
231 : output_format=outformat, &
232 : content=content_id(id), &
233 : title=title(id), &
234 : cell=cell, &
235 : unit_conv=unit_conv, &
236 468 : print_kind=print_kind)
237 :
238 : CALL cp_print_key_finished_output(unit_nr, logger, &
239 468 : print_key, "", local=.FALSE.)
240 :
241 : END IF
242 :
243 964 : CALL f_env_rm_defaults(f_env, ierr, handle)
244 3374 : CPASSERT(ierr == 0)
245 :
246 : END DO
247 :
248 482 : CALL timestop(handle1)
249 482 : END SUBROUTINE pint_write_centroids
250 :
251 : ! ***************************************************************************
252 : !> \brief Write out the trajectory of the beads (positions and velocities)
253 : !> \param pint_env ...
254 : !> \par History
255 : !> 2010-11-25 added support for velocity printing [lwalewski]
256 : !> \author hforbert
257 : ! **************************************************************************************************
258 482 : SUBROUTINE pint_write_trajectory(pint_env)
259 : TYPE(pint_env_type), INTENT(IN) :: pint_env
260 :
261 : CHARACTER(len=*), PARAMETER :: routineN = 'pint_write_trajectory'
262 : INTEGER, PARAMETER :: force_id = 3, n_ids = 3, pos_id = 1, &
263 : vel_id = 2
264 :
265 : CHARACTER(len=default_string_length) :: ext, form, ib_str, my_middle_name, &
266 : title, unit_str
267 : CHARACTER(len=default_string_length), DIMENSION(3) :: content_id, middle_name, sect_path
268 : INTEGER :: handle, handle1, iat, ib, id, idim, &
269 : idir, ierr, imag_stride, outformat, &
270 : should_output, unit_nr
271 : LOGICAL :: new_file
272 : REAL(kind=dp) :: unit_conv
273 : TYPE(cell_type), POINTER :: cell
274 : TYPE(cp_logger_type), POINTER :: logger
275 : TYPE(cp_subsys_type), POINTER :: subsys
276 : TYPE(f_env_type), POINTER :: f_env
277 : TYPE(particle_list_type), POINTER :: particles
278 : TYPE(section_vals_type), POINTER :: print_key
279 :
280 482 : CALL timeset(routineN, handle1)
281 :
282 482 : sect_path(pos_id) = "MOTION%PRINT%TRAJECTORY"
283 482 : sect_path(vel_id) = "MOTION%PRINT%VELOCITIES"
284 482 : sect_path(force_id) = "MOTION%PRINT%FORCES"
285 482 : middle_name(pos_id) = "pos-"
286 482 : middle_name(vel_id) = "vel-"
287 482 : middle_name(force_id) = "force-"
288 482 : content_id(pos_id) = "POS"
289 482 : content_id(vel_id) = "VEL"
290 482 : content_id(force_id) = "FORCE"
291 :
292 482 : NULLIFY (logger)
293 482 : logger => cp_get_default_logger()
294 :
295 482 : CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
296 :
297 : ! iterate over the properties that we know how to print
298 : ! (currently positions and velocities)
299 1928 : DO id = 1, n_ids
300 :
301 : print_key => section_vals_get_subs_vals(pint_env%input, &
302 1446 : TRIM(sect_path(id)))
303 :
304 : should_output = cp_print_key_should_output( &
305 : iteration_info=logger%iter_info, &
306 1446 : basis_section=print_key)
307 : IF (.NOT. BTEST(should_output, cp_p_file)) CONTINUE
308 :
309 : ! get units of measure for output (if available)
310 : CALL section_vals_val_get(print_key, "UNIT", &
311 1446 : c_val=unit_str)
312 1446 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
313 :
314 : ! get the format for output
315 1446 : CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
316 :
317 0 : SELECT CASE (outformat)
318 : CASE (dump_dcd, dump_dcd_aligned_cell)
319 0 : form = "UNFORMATTED"
320 0 : ext = ".dcd"
321 : CASE (dump_atomic)
322 0 : form = "FORMATTED"
323 0 : ext = ""
324 : CASE (dump_xmol)
325 1446 : form = "FORMATTED"
326 1446 : ext = ".xyz"
327 : CASE default
328 1446 : CPABORT("")
329 : END SELECT
330 :
331 1446 : NULLIFY (f_env, cell, subsys)
332 : CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
333 1446 : f_env=f_env, handle=handle)
334 : CALL force_env_get(force_env=f_env%force_env, &
335 1446 : cell=cell, subsys=subsys)
336 1446 : CALL cp_subsys_get(subsys, particles=particles)
337 :
338 : !Get print stride for bead trajectories
339 : CALL section_vals_val_get(pint_env%input, &
340 : "MOTION%PINT%PRINT%IMAGINARY_TIME_STRIDE", &
341 1446 : i_val=imag_stride)
342 :
343 : ! iterate over beads
344 11334 : DO ib = 1, pint_env%p, imag_stride
345 :
346 : ! copy the requested property of the current bead
347 : ! to the particles structure
348 9888 : idim = 0
349 1696776 : DO iat = 1, pint_env%ndim/3
350 6757440 : DO idir = 1, 3
351 5060664 : idim = idim + 1
352 5060664 : particles%els(iat)%r(idir) = pint_env%x(ib, idim)
353 5060664 : particles%els(iat)%v(idir) = pint_env%v(ib, idim)
354 6747552 : particles%els(iat)%f(idir) = pint_env%f(ib, idim)
355 : END DO
356 : END DO
357 :
358 : ! set up the output unit number and file name
359 : ! for the current property and bead
360 9888 : ib_str = ""
361 9888 : WRITE (ib_str, *) ib
362 9888 : my_middle_name = TRIM(middle_name(id))//TRIM(ADJUSTL(ib_str))
363 : unit_nr = cp_print_key_unit_nr(logger=logger, &
364 : basis_section=print_key, print_key_path="", &
365 : extension=TRIM(ext), middle_name=TRIM(my_middle_name), &
366 9888 : local=.FALSE., file_form=form, is_new_file=new_file)
367 :
368 : ! don't write the 0-th frame if the file already exists
369 9888 : IF (.NOT. new_file .AND. (pint_env%iter .LE. pint_env%first_step)) THEN
370 : CALL cp_print_key_finished_output(unit_nr, logger, &
371 808 : print_key)
372 : CONTINUE
373 : END IF
374 :
375 : ! actually perform the i/o - on the ionode only
376 11334 : IF (unit_nr > 0) THEN
377 :
378 1620 : IF (outformat == dump_xmol) THEN
379 : WRITE (UNIT=title, FMT="(A,I8,A,F20.10)") &
380 1620 : " i =", pint_env%iter, &
381 3240 : ", E =", pint_env%e_pot_bead(ib)
382 : END IF
383 :
384 : CALL write_particle_coordinates( &
385 : particles%els, &
386 : iunit=unit_nr, &
387 : output_format=outformat, &
388 : content=content_id(id), &
389 : title=title, &
390 : cell=cell, &
391 1620 : unit_conv=unit_conv)
392 :
393 : CALL cp_print_key_finished_output(unit_nr, logger, &
394 1620 : print_key, "", local=.FALSE.)
395 :
396 : END IF
397 :
398 : END DO
399 :
400 1446 : CALL f_env_rm_defaults(f_env, ierr, handle)
401 6266 : CPASSERT(ierr == 0)
402 :
403 : END DO
404 :
405 482 : CALL timestop(handle1)
406 482 : END SUBROUTINE pint_write_trajectory
407 :
408 : ! ***************************************************************************
409 : !> \brief Write center of mass (COM) position according to PINT%PRINT%COM
410 : !> \param pint_env ...
411 : !> \date 2010-02-17
412 : !> \author Lukasz Walewski
413 : ! **************************************************************************************************
414 482 : SUBROUTINE pint_write_com(pint_env)
415 :
416 : TYPE(pint_env_type), INTENT(IN) :: pint_env
417 :
418 : CHARACTER(len=default_string_length) :: stmp1, stmp2
419 : INTEGER :: ic, unit_nr
420 : LOGICAL :: new_file, should_output
421 : REAL(kind=dp), DIMENSION(3) :: com_r
422 : TYPE(cp_logger_type), POINTER :: logger
423 : TYPE(section_vals_type), POINTER :: print_key
424 :
425 482 : NULLIFY (logger)
426 482 : logger => cp_get_default_logger()
427 :
428 : ! decide whether to write anything or not
429 482 : NULLIFY (print_key)
430 : print_key => section_vals_get_subs_vals(pint_env%input, &
431 482 : "MOTION%PINT%PRINT%COM")
432 : should_output = BTEST(cp_print_key_should_output( &
433 : iteration_info=logger%iter_info, &
434 482 : basis_section=print_key), cp_p_file)
435 482 : IF (.NOT. should_output) THEN
436 482 : RETURN
437 : END IF
438 :
439 0 : com_r = pint_com_pos(pint_env)
440 0 : DO ic = 1, 3
441 0 : com_r(ic) = cp_unit_from_cp2k(com_r(ic), "angstrom")
442 : END DO
443 :
444 : unit_nr = cp_print_key_unit_nr(logger, print_key, is_new_file=new_file, &
445 0 : middle_name="com-pos", extension=".xyz")
446 :
447 : ! don't write the 0-th frame if the file already exists
448 0 : IF (.NOT. new_file .AND. (pint_env%iter .LE. pint_env%first_step)) THEN
449 : CALL cp_print_key_finished_output(unit_nr, logger, &
450 0 : print_key)
451 0 : RETURN
452 : END IF
453 :
454 : ! actually perform the i/o - on the ionode only
455 0 : IF (unit_nr > 0) THEN
456 :
457 0 : WRITE (unit_nr, '(I2)') 1
458 0 : WRITE (stmp1, *) pint_env%iter
459 0 : WRITE (stmp2, '(F20.10)') pint_env%energy(e_conserved_id)
460 0 : WRITE (unit_nr, '(4A)') " Iteration = ", TRIM(ADJUSTL(stmp1)), &
461 0 : ", E_conserved = ", TRIM(ADJUSTL(stmp2))
462 0 : WRITE (unit_nr, '(A2,3(1X,F20.10))') "X ", (com_r(ic), ic=1, 3)
463 :
464 0 : CALL m_flush(unit_nr)
465 :
466 : END IF
467 :
468 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
469 :
470 : END SUBROUTINE pint_write_com
471 :
472 : ! ***************************************************************************
473 : !> \brief Writes out the energies according to PINT%PRINT%ENERGY
474 : !> \param pint_env path integral environment
475 : !> \par History
476 : !> various bug fixes [hforbert]
477 : !> 2009-11-16 energy components calc moved out of here [lwalewski]
478 : !> \author fawzi
479 : ! **************************************************************************************************
480 482 : SUBROUTINE pint_write_ener(pint_env)
481 : TYPE(pint_env_type), INTENT(IN) :: pint_env
482 :
483 : INTEGER :: ndof, unit_nr
484 : LOGICAL :: file_is_new
485 : REAL(kind=dp) :: t, temp
486 : TYPE(cp_logger_type), POINTER :: logger
487 : TYPE(section_vals_type), POINTER :: print_key
488 :
489 482 : NULLIFY (print_key, logger)
490 : print_key => section_vals_get_subs_vals(pint_env%input, &
491 482 : "MOTION%PINT%PRINT%ENERGY")
492 482 : logger => cp_get_default_logger()
493 482 : IF (BTEST(cp_print_key_should_output(iteration_info=logger%iter_info, &
494 : basis_section=print_key), cp_p_file)) THEN
495 :
496 : unit_nr = cp_print_key_unit_nr(logger, print_key, middle_name="energy", &
497 482 : extension=".dat", is_new_file=file_is_new)
498 :
499 : ! don't write the 0-th frame if the file already exists
500 482 : IF (.NOT. file_is_new .AND. (pint_env%iter .LE. pint_env%first_step)) THEN
501 : CALL cp_print_key_finished_output(unit_nr, logger, &
502 34 : print_key)
503 34 : RETURN
504 : END IF
505 :
506 : ! cp_print_key_unit_nr returns -1 on nodes other than logger%para_env%is_source()
507 448 : IF (unit_nr > 0) THEN
508 :
509 : ! please keep the format explanation up to date
510 : ! keep the constant of motion the true constant of motion !
511 234 : IF (file_is_new) THEN
512 : WRITE (unit_nr, "(A8,1X,A12,1X,5(A20,1X),A12)") &
513 20 : "# StepNr", &
514 20 : " Time [fs]", &
515 20 : " Kinetic [a.u.]", &
516 20 : " VirialKin [a.u.]", &
517 20 : " Temperature [K]", &
518 20 : " Potential [a.u.]", &
519 20 : " ConsQty [a.u.]", &
520 40 : " CPU [s]"
521 : END IF
522 :
523 234 : t = cp_unit_from_cp2k(pint_env%t, "fs")
524 :
525 234 : ndof = pint_env%p
526 234 : IF (pint_env%first_propagated_mode .EQ. 2) THEN
527 0 : ndof = ndof - 1
528 : END IF
529 : temp = cp_unit_from_cp2k(2.0_dp*pint_env%e_kin_beads/ &
530 : REAL(ndof, dp)/REAL(pint_env%ndim, dp), &
531 234 : "K")*pint_env%propagator%temp_sim2phys
532 :
533 : WRITE (unit_nr, "(I8,1X,F12.3,1X,5(F20.9,1X),F12.1)") &
534 234 : pint_env%iter, &
535 234 : t, &
536 234 : pint_env%energy(e_kin_thermo_id), &
537 234 : pint_env%energy(e_kin_virial_id), &
538 234 : temp, &
539 234 : pint_env%energy(e_potential_id), &
540 234 : pint_env%energy(e_conserved_id), &
541 468 : pint_env%time_per_step
542 234 : CALL m_flush(unit_nr)
543 :
544 : END IF
545 :
546 448 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
547 : END IF
548 :
549 : END SUBROUTINE pint_write_ener
550 :
551 : ! ***************************************************************************
552 : !> \brief Writes out the actions according to PINT%PRINT%ACTION
553 : !> \param pint_env path integral environment
554 : !> \author Felix Uhl
555 : ! **************************************************************************************************
556 482 : SUBROUTINE pint_write_action(pint_env)
557 : TYPE(pint_env_type), INTENT(IN) :: pint_env
558 :
559 : INTEGER :: unit_nr
560 : LOGICAL :: file_is_new
561 : REAL(kind=dp) :: t
562 : TYPE(cp_logger_type), POINTER :: logger
563 : TYPE(section_vals_type), POINTER :: print_key
564 :
565 482 : NULLIFY (print_key, logger)
566 : print_key => section_vals_get_subs_vals(pint_env%input, &
567 482 : "MOTION%PINT%PRINT%ACTION")
568 482 : logger => cp_get_default_logger()
569 482 : IF (BTEST(cp_print_key_should_output(iteration_info=logger%iter_info, &
570 : basis_section=print_key), cp_p_file)) THEN
571 :
572 : unit_nr = cp_print_key_unit_nr(logger, print_key, middle_name="action", &
573 148 : extension=".dat", is_new_file=file_is_new)
574 :
575 : ! don't write the 0-th frame if the file already exists
576 148 : IF (.NOT. file_is_new .AND. (pint_env%iter .LE. pint_env%first_step)) THEN
577 : CALL cp_print_key_finished_output(unit_nr, logger, &
578 16 : print_key)
579 16 : RETURN
580 : END IF
581 :
582 : ! cp_print_key_unit_nr returns -1 on nodes other than logger%para_env%is_source()
583 132 : IF (unit_nr > 0) THEN
584 :
585 : ! please keep the format explanation up to date
586 : ! keep the constant of motion the true constant of motion !
587 72 : IF (file_is_new) THEN
588 : WRITE (unit_nr, "(A8,1X,A12,1X,2(A25,1X))") &
589 12 : "# StepNr", &
590 12 : " Time [fs]", &
591 12 : " Link Action [a.u.]", &
592 24 : " Potential Action [a.u.]"
593 : END IF
594 :
595 72 : t = cp_unit_from_cp2k(pint_env%t, "fs")
596 :
597 : WRITE (unit_nr, "(I8,1X,F12.3,1X,5(F20.9,1X),F12.1)") &
598 72 : pint_env%iter, &
599 72 : t, &
600 72 : pint_env%link_action, &
601 144 : pint_env%pot_action
602 72 : CALL m_flush(unit_nr)
603 :
604 : END IF
605 :
606 132 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
607 : END IF
608 :
609 : END SUBROUTINE pint_write_action
610 :
611 : ! ***************************************************************************
612 : !> \brief Write step info to the output file.
613 : !> \param pint_env ...
614 : !> \date 2009-11-16
615 : !> \par History
616 : !> 2010-01-27 getting default unit nr now only on ionode [lwalewski]
617 : !> \author Lukasz Walewski
618 : ! **************************************************************************************************
619 428 : SUBROUTINE pint_write_step_info(pint_env)
620 : TYPE(pint_env_type), INTENT(IN) :: pint_env
621 :
622 : CHARACTER(len=default_string_length) :: msgstr, stmp, time_unit
623 : INTEGER :: unit_nr
624 : REAL(kind=dp) :: time_used
625 : TYPE(cp_logger_type), POINTER :: logger
626 :
627 428 : unit_nr = 0
628 428 : NULLIFY (logger)
629 428 : logger => cp_get_default_logger()
630 :
631 428 : time_used = pint_env%time_per_step
632 428 : time_unit = "sec"
633 428 : IF (time_used .GE. 60.0_dp) THEN
634 0 : time_used = time_used/60.0_dp
635 0 : time_unit = "min"
636 : END IF
637 428 : IF (time_used .GE. 60.0_dp) THEN
638 0 : time_used = time_used/60.0_dp
639 0 : time_unit = "hours"
640 : END IF
641 428 : msgstr = "PINT step"
642 428 : stmp = ""
643 428 : WRITE (stmp, *) pint_env%iter
644 428 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))//" of"
645 428 : stmp = ""
646 428 : WRITE (stmp, *) pint_env%last_step
647 428 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))//" in"
648 428 : stmp = ""
649 428 : WRITE (stmp, '(F20.1)') time_used
650 428 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))
651 428 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(time_unit))//"."
652 :
653 428 : IF (logger%para_env%is_source()) THEN
654 214 : unit_nr = cp_logger_get_default_unit_nr(logger)
655 214 : WRITE (unit_nr, '(T2,A)') "PINT| "//TRIM(ADJUSTL(msgstr))
656 : END IF
657 :
658 : ! print out the total energy - for regtest evaluation
659 428 : stmp = ""
660 428 : WRITE (stmp, *) pint_env%energy(e_conserved_id)
661 428 : msgstr = "Total energy = "//TRIM(ADJUSTL(stmp))
662 428 : IF (logger%para_env%is_source()) THEN
663 214 : WRITE (unit_nr, '(T2,A)') "PINT| "//TRIM(ADJUSTL(msgstr))
664 : END IF
665 :
666 428 : END SUBROUTINE pint_write_step_info
667 :
668 : ! ***************************************************************************
669 : !> \brief Write radii of gyration according to PINT%PRINT%CENTROID_GYR
670 : !> \param pint_env ...
671 : !> \date 2011-01-07
672 : !> \author Lukasz Walewski
673 : ! **************************************************************************************************
674 482 : SUBROUTINE pint_write_rgyr(pint_env)
675 :
676 : TYPE(pint_env_type), INTENT(IN) :: pint_env
677 :
678 : CHARACTER(len=default_string_length) :: unit_str
679 : INTEGER :: ia, ib, ic, idim, unit_nr
680 : LOGICAL :: new_file, should_output
681 : REAL(kind=dp) :: nb, ss, unit_conv
682 : TYPE(cp_logger_type), POINTER :: logger
683 : TYPE(section_vals_type), POINTER :: print_key
684 :
685 482 : NULLIFY (logger)
686 482 : logger => cp_get_default_logger()
687 :
688 : ! decide whether to write anything or not
689 482 : NULLIFY (print_key)
690 : print_key => section_vals_get_subs_vals(pint_env%input, &
691 482 : "MOTION%PINT%PRINT%CENTROID_GYR")
692 : should_output = BTEST(cp_print_key_should_output( &
693 : iteration_info=logger%iter_info, &
694 482 : basis_section=print_key), cp_p_file)
695 482 : IF (.NOT. should_output) THEN
696 98 : RETURN
697 : END IF
698 :
699 : ! get the units conversion factor
700 410 : CALL section_vals_val_get(print_key, "UNIT", c_val=unit_str)
701 410 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
702 :
703 : ! calculate the centroid positions
704 410 : nb = REAL(pint_env%p, dp)
705 410 : idim = 0
706 1706 : DO ia = 1, pint_env%ndim/3
707 5594 : DO ic = 1, 3
708 3888 : idim = idim + 1
709 3888 : ss = 0.0_dp
710 31896 : DO ib = 1, pint_env%p
711 31896 : ss = ss + pint_env%x(ib, idim)
712 : END DO
713 5184 : pint_env%rtmp_ndim(idim) = ss/nb
714 : END DO
715 : END DO
716 :
717 : ! calculate the radii of gyration
718 : idim = 0
719 1706 : DO ia = 1, pint_env%ndim/3
720 : ss = 0.0_dp
721 5184 : DO ic = 1, 3
722 3888 : idim = idim + 1
723 33192 : DO ib = 1, pint_env%p
724 31896 : ss = ss + (pint_env%x(ib, idim) - pint_env%rtmp_ndim(idim))**2
725 : END DO
726 : END DO
727 1706 : pint_env%rtmp_natom(ia) = SQRT(ss/nb)*unit_conv
728 : END DO
729 :
730 : unit_nr = cp_print_key_unit_nr(logger, print_key, is_new_file=new_file, &
731 410 : middle_name="centroid-gyr", extension=".dat")
732 :
733 : ! don't write the 0-th frame if the file already exists
734 410 : IF (.NOT. new_file .AND. (pint_env%iter .LE. pint_env%first_step)) THEN
735 : CALL cp_print_key_finished_output(unit_nr, logger, &
736 26 : print_key)
737 26 : RETURN
738 : END IF
739 :
740 : ! actually perform the i/o - on the ionode only
741 384 : IF (unit_nr > 0) THEN
742 :
743 829 : DO ia = 1, pint_env%ndim/3
744 829 : WRITE (unit_nr, '(F20.10,1X)', ADVANCE='NO') pint_env%rtmp_natom(ia)
745 : END DO
746 199 : WRITE (unit_nr, '(A)') ""
747 :
748 199 : CALL m_flush(unit_nr)
749 :
750 : END IF
751 :
752 384 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
753 :
754 : END SUBROUTINE pint_write_rgyr
755 :
756 : END MODULE pint_io
|
