(git:98357aa)
Loading...
Searching...
No Matches
dumpdcd.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8PROGRAM dumpdcd
9
10! Version: 3.2
11! Author: Matthias Krack (MK)
12! History: - Creation (13.02.2012,MK)
13! - XYZ file option added (13.02.2012,MK)
14! - Flags added for first and last frame and -o for output (24.05.2012,MK)
15! - XMOL flag added (25.05.2012,MK)
16! - PBC flag added (04.06.2012,MK)
17! - Stride flag added (05.06.2012,MK)
18! - Tracing of atoms added to detect the atoms which left the box (06.06.2012,MK)
19! - Keep input coordinates for further processing steps (15.06.2012,MK)
20! - VEL to CORD (-vel2cord flag) hack added (25.06.2012,MK)
21! - Added -displacement (-disp) flag (26.06.2012,MK)
22! - Dump the atomic displacements (CORD file) or temperatures (VEL file as x-coordinates of a DCD file (28.06.2012,MK)
23! - FRC to CORD (-frc2cord flag) hack added (28.01.2016,MK)
24! - Option -eformat added (15.02.2016,MK)
25! - Fix box unit string for velocity dump (07.06.2020,MK)
26! - Dump cell information to comment line of XMOL file if requested and available (15.02.2021,MK)
27! - Write comment line compliant with REFTRAJ format (16.06.2022,MK)
28!
29! Note: For -ekin a XYZ file is required to obtain the atomic labels.
30! The -info and the -debug flags provide a more detailed output which is especially handy for tracing problems.
31! The output in DCD format is in binary format.
32! The input coordinates should be in Angstrom. Velocities and forces are expected to be in atomic units.
33
34! Uncomment the following line if this module is available (e.g. with gfortran) and comment the corresponding
35! variable declarations below
36! USE ISO_FORTRAN_ENV, ONLY: error_unit,input_unit,output_unit
37
38 IMPLICIT NONE
39
40 ! Comment the following lines if the ISO_FORTRAN_ENV is used (see above)
41 INTEGER, PARAMETER :: default_error_unit = 0, &
42 default_input_unit = 5, &
43 default_output_unit = 6
44 INTEGER :: error_unit = default_error_unit, &
45 input_unit = default_input_unit, &
46 output_unit = default_output_unit
47 ! End Comment
48
49 ! Parameters
50 CHARACTER(LEN=*), PARAMETER :: routinen = "dumpdcd", &
51 version_info = routinen//" v3.2 (16.06.2022)"
52
53 INTEGER, PARAMETER :: dp = selected_real_kind(14, 200), &
54 sp = selected_real_kind(6, 30)
55 INTEGER, PARAMETER :: default_string_length = 240, &
56 cell_input_unit = 10, &
57 xyz_input_unit = 11
58
59 REAL(kind=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
60 REAL(kind=dp), PARAMETER :: angstrom = 0.52917720859_dp, & ! [a.u.] -> [Angstrom]
61 degree = 180.0_dp/pi, & ! [rad] -> [degree]
62 kelvin = 315774.647902944_dp, & ! [a.u.] -> [K]
63 massunit = 1822.88484264550_dp ! [u] -> [a.u.]
64
65 ! Variables
66 CHARACTER(LEN=4) :: id_dcd
67 CHARACTER(LEN=10) :: unit_string
68 CHARACTER(LEN=17) :: fmt_string
69 CHARACTER(LEN=default_string_length) :: arg, cell_file_name, dcd_file_name, message, out_file_name, &
70 output_format, remark_xyz, string, xyz_file_name
71 CHARACTER(LEN=5), DIMENSION(:), ALLOCATABLE :: atomic_label
72 CHARACTER(LEN=80), DIMENSION(2) :: remark_dcd
73 INTEGER :: first_frame, have_unit_cell, i, iarg, iatom, &
74 iframe, istat, istep_dcd, istride_dcd, last_frame, narg, &
75 natom_dcd, natom_xyz, ndcd_file, nframe, &
76 nframe_read, nremark, stride
77 LOGICAL :: apply_pbc, debug, dump_frame, eformat, ekin, eo, &
78 have_atomic_labels, have_cell_file, ignore_warnings, info, &
79 opened, output_format_dcd, output_format_xmol, pbc0, &
80 print_atomic_displacements, print_scaled_coordinates, &
81 print_scaled_pbc_coordinates, trace_atoms, vel2cord
82 REAL(kind=sp) :: dt_dcd
83 REAL(kind=dp) :: a, a_dcd, alpha, alpha_dcd, b, b_dcd, beta, beta_dcd, c, c_dcd, &
84 cell_volume, dt, energy, eps_angle, eps_geo, eps_out_of_box, &
85 first_step_time, gamma, gamma_dcd, md_time_step, ndt, step_time, &
86 tavg, tavg_frame
87 INTEGER, DIMENSION(16) :: idum
88 REAL(kind=dp), DIMENSION(3) :: rdum
89 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: atomic_displacement, atomic_mass, atomic_temperature
90 REAL(kind=dp), DIMENSION(3, 3) :: h, hinv, hmat
91 REAL(kind=sp), DIMENSION(:, :), ALLOCATABLE :: r
92 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: r_pbc, r0, s, s_pbc
93
94 apply_pbc = .false.
95 pbc0 = .false.
96 debug = .false.
97 dump_frame = .true.
98 eformat = .false.
99 ekin = .false.
100 eo = .false.
101 ignore_warnings = .false.
102 info = .false.
103 trace_atoms = .false.
104 vel2cord = .false.
105 print_atomic_displacements = .false.
106 print_scaled_coordinates = .false.
107 print_scaled_pbc_coordinates = .false.
108 first_frame = 1
109 last_frame = 1000000 ! Hard limit of 1 Mio frames in total
110 stride = 1
111 ndcd_file = 0
112 nframe = 0
113 nframe_read = 0
114 nremark = 0
115 have_unit_cell = 0
116 have_atomic_labels = .false.
117 have_cell_file = .false.
118 idum(:) = 0
119 rdum(:) = 0.0_dp
120 id_dcd = ""
121 cell_file_name = ""
122 dcd_file_name = ""
123 out_file_name = ""
124 xyz_file_name = ""
125 output_format = "default"
126 output_format_dcd = .false.
127 output_format_xmol = .false.
128 remark_dcd(:) = ""
129 remark_xyz = ""
130 fmt_string = ""
131 dt = 0.0_dp
132 dt_dcd = 0.0_sp
133 a = 0.0_dp
134 b = 0.0_dp
135 c = 0.0_dp
136 alpha = 0.0_dp
137 beta = 0.0_dp
138 gamma = 0.0_dp
139 eps_out_of_box = -huge(0.0_dp)
140 first_step_time = 0.0_dp
141 md_time_step = 0.0_dp
142 step_time = 0.0_dp
143 tavg = 0.0_dp
144 tavg_frame = 0.0_dp
145 energy = 0.0_dp
146
147 narg = command_argument_count()
148
149 IF (narg == 0) THEN
150 CALL print_help()
151 CALL abort_program(routinen, "No input file(s) specified")
152 END IF
153
154 iarg = 0
155
156 dcd_file_loop: DO
157
158 iarg = iarg + 1
159
160 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
161
162 SELECT CASE (arg)
163 CASE ("-cell_file", "-cell")
164 iarg = iarg + 1
165 CALL get_command_argument(number=iarg, VALUE=cell_file_name, status=istat)
166 have_cell_file = .true.
167 cycle dcd_file_loop
168 CASE ("-debug", "-d")
169 debug = .true.
170 info = .true.
171 cycle dcd_file_loop
172 CASE ("-displacements", "-disp")
173 print_atomic_displacements = .true.
174 cycle dcd_file_loop
175 CASE ("-eformat")
176 eformat = .true.
177 cycle dcd_file_loop
178 CASE ("-ekin")
179 ekin = .true.
180 cycle dcd_file_loop
181 CASE ("-eo")
182 eo = .true.
183 cycle dcd_file_loop
184 CASE ("-first_frame", "-first", "-ff")
185 iarg = iarg + 1
186 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
187 READ (unit=arg, fmt=*, err=100) first_frame
188 IF (first_frame <= 0) THEN
189 CALL abort_program(routinen, "Invalid number for first frame specified: "// &
190 "first_frame must be greater than zero")
191 END IF
192 cycle dcd_file_loop
193100 CALL abort_program(routinen, "Invalid number for first frame specified "// &
194 "(an integer number greater than zero is expected)")
195 CASE ("-help", "-h")
196 CALL print_help()
197 stop
198 CASE ("-ignore_warnings")
199 ignore_warnings = .true.
200 cycle dcd_file_loop
201 CASE ("-info", "-i")
202 info = .true.
203 cycle dcd_file_loop
204 CASE ("-last_frame", "-last", "-lf")
205 iarg = iarg + 1
206 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
207 READ (unit=arg, fmt=*, err=101) last_frame
208 IF (last_frame <= 0) THEN
209 CALL abort_program(routinen, "Invalid number for last frame specified: "// &
210 "last_frame must be greater than zero")
211 END IF
212 cycle dcd_file_loop
213101 CALL abort_program(routinen, "Invalid number for last frame specified "// &
214 "(an integer number greater than zero is expected)")
215 CASE ("-md_time_step")
216 iarg = iarg + 1
217 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
218 READ (unit=arg, fmt=*, err=102) md_time_step
219 IF (md_time_step <= 0.0_dp) THEN
220 CALL abort_program(routinen, "Invalid (negative) MD time step specified")
221 END IF
222 cycle dcd_file_loop
223102 CALL abort_program(routinen, "Invalid MD time step specified")
224 CASE ("-o", "-output")
225 iarg = iarg + 1
226 CALL get_command_argument(number=iarg, VALUE=out_file_name, status=istat)
227 cycle dcd_file_loop
228 CASE ("-output_format", "-of")
229 iarg = iarg + 1
230 CALL get_command_argument(number=iarg, VALUE=output_format, status=istat)
231 CALL uppercase(output_format)
232 SELECT CASE (output_format)
233 CASE ("DCD")
234 output_format_dcd = .true.
235 output_format_xmol = .false.
236 CASE ("XMOL")
237 output_format_dcd = .false.
238 output_format_xmol = .true.
239 CASE DEFAULT
240 CALL abort_program(routinen, "Invalid output format type specified")
241 END SELECT
242 cycle dcd_file_loop
243 CASE ("-pbc")
244 apply_pbc = .true.
245 pbc0 = .false.
246 cycle dcd_file_loop
247 CASE ("-pbc0")
248 apply_pbc = .true.
249 pbc0 = .true.
250 cycle dcd_file_loop
251 CASE ("-scaled_coordinates", "-sc")
252 print_scaled_coordinates = .true.
253 cycle dcd_file_loop
254 CASE ("-scaled_pbc_coordinates", "-spc")
255 print_scaled_pbc_coordinates = .true.
256 cycle dcd_file_loop
257 CASE ("-stride")
258 iarg = iarg + 1
259 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
260 READ (unit=arg, fmt=*, err=104) stride
261 IF (stride < 1) THEN
262 CALL abort_program(routinen, "Invalid stride for frame dump specified: stride must be greater than zero")
263 END IF
264 cycle dcd_file_loop
265104 CALL abort_program(routinen, "Invalid stride for frame dump specified "// &
266 "(an integer number greater than 0 is expected)")
267 CASE ("-trace_atoms")
268 iarg = iarg + 1
269 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
270 READ (unit=arg, fmt=*, err=108) eps_out_of_box
271 IF (eps_out_of_box <= 0.0_dp) THEN
272 CALL abort_program(routinen, "Invalid threshold value for -trace_atoms flag specified")
273 END IF
274 trace_atoms = .true.
275 cycle dcd_file_loop
276108 CALL abort_program(routinen, "Invalid threshold value for -trace_atoms flag specified")
277 CASE ("-vel2cord", "-v2c", "-frc2cord", "-f2c")
278 vel2cord = .true.
279 cycle dcd_file_loop
280 CASE ("-xyz", "xyz_file")
281 iarg = iarg + 1
282 CALL get_command_argument(number=iarg, VALUE=xyz_file_name, status=istat)
283 have_atomic_labels = .true.
284 cycle dcd_file_loop
285 CASE DEFAULT
286 IF (arg(1:1) == "-") THEN
287 CALL print_help()
288 CALL abort_program(routinen, "Unknown command line flag """//trim(arg)//""" found")
289 END IF
290 dcd_file_name = arg
291 END SELECT
292
293 ! Check flag compatibility
294 IF (first_frame > last_frame) THEN
295 CALL abort_program(routinen, "Number of first frame greater than number of last frame")
296 END IF
297 IF ((.NOT. have_atomic_labels) .AND. output_format_xmol) THEN
298 CALL abort_program(routinen, "The output format XMOL requires a valid xyz file (-xyz flag)")
299 END IF
300 IF (output_format_xmol .AND. ekin) THEN
301 CALL abort_program(routinen, "Output format XMOL and the -ekin flag are incompatible")
302 END IF
303 IF (output_format_xmol .AND. print_atomic_displacements) THEN
304 CALL abort_program(routinen, "Output format XMOL and the -displacements flag are incompatible")
305 END IF
306 IF (output_format_xmol .AND. eo) THEN
307 CALL abort_program(routinen, "The -eo flag is incompatible with the output format XMOL")
308 END IF
309 IF (.NOT. have_atomic_labels .AND. ekin) THEN
310 CALL abort_program(routinen, "ekin flag requires also the specification of a valid xyz file (-xyz flag)")
311 END IF
312 IF (.NOT. apply_pbc .AND. trace_atoms) THEN
313 CALL abort_program(routinen, "The -trace_atoms flag requires the specification of a -pbc flag")
314 END IF
315 IF (ekin .AND. print_atomic_displacements) THEN
316 CALL abort_program(routinen, "The -ekin flag and the -displacements flag are incompatible")
317 END IF
318 IF (print_scaled_coordinates .AND. print_scaled_pbc_coordinates) THEN
319 CALL abort_program(routinen, "The -sc flag and the -spc flag are incompatible")
320 END IF
321 IF (.NOT. apply_pbc .AND. print_scaled_coordinates) THEN
322 CALL abort_program(routinen, "The -sc flag requires the specification of a -pbc flag")
323 END IF
324 IF (.NOT. apply_pbc .AND. print_scaled_pbc_coordinates) THEN
325 CALL abort_program(routinen, "The -spc flag requires the specification of a -pbc flag")
326 END IF
327
328 ! Use optionally an ES format
329 IF (eformat) THEN
330 fmt_string = "(A5,3(1X,ES14.6))"
331 ELSE
332 fmt_string = "(A5,3(1X, F14.6))"
333 END IF
334
335 ! Open output units as requested
336 IF (output_format_dcd) THEN
337 ! Set default output file name if no file name was specified
338 IF (len_trim(out_file_name) == 0) out_file_name = "output.dcd"
339 ! Check if a new output file name was specified, if yes then close the old unit
340 INQUIRE (unit=output_unit, name=string)
341 IF (trim(string) /= trim(out_file_name)) CLOSE (unit=output_unit)
342 INQUIRE (unit=output_unit, opened=opened)
343 IF (.NOT. opened) THEN
344 OPEN (unit=output_unit, &
345 file=out_file_name, &
346 status="UNKNOWN", &
347 access="SEQUENTIAL", &
348 form="UNFORMATTED", &
349 action="WRITE", &
350 iostat=istat)
351 IF (istat /= 0) CALL abort_program(routinen, "The unformatted output file could not be opened")
352 END IF
353 ELSE
354 IF (eo) error_unit = output_unit
355 IF (len_trim(out_file_name) > 0) THEN
356 ! Check if a new output file name was specified, if yes then close the old unit
357 INQUIRE (unit=output_unit, name=string)
358 IF (trim(string) /= trim(out_file_name)) CLOSE (unit=output_unit)
359 INQUIRE (unit=output_unit, opened=opened)
360 IF (.NOT. opened) THEN
361 OPEN (unit=output_unit, &
362 file=out_file_name, &
363 status="UNKNOWN", &
364 access="SEQUENTIAL", &
365 form="FORMATTED", &
366 position="REWIND", &
367 action="WRITE", &
368 iostat=istat)
369 IF (istat /= 0) CALL abort_program(routinen, "The formatted output file could not be opened")
370 END IF
371 END IF
372 END IF
373
374 ! Avoid reading from and writing to the same file
375 IF (trim(dcd_file_name) == trim(out_file_name)) THEN
376 CALL abort_program(routinen, "Input and output file name cannot be the same")
377 END IF
378
379 ! Read reference xyz input file if requested
380 IF (have_atomic_labels) THEN
381 string = ""
382 ! Check if a new XYZ input file name was specified, if yes then close the old unit
383 INQUIRE (unit=xyz_input_unit, name=string)
384 IF (trim(string) /= trim(xyz_file_name)) CLOSE (unit=xyz_input_unit)
385 INQUIRE (unit=xyz_input_unit, opened=opened)
386 ! Read new XYZ input file if needed and update reference information
387 IF (.NOT. opened) THEN
388 OPEN (unit=xyz_input_unit, &
389 file=xyz_file_name, &
390 status="OLD", &
391 access="SEQUENTIAL", &
392 form="FORMATTED", &
393 position="REWIND", &
394 action="READ", &
395 iostat=istat)
396 IF (istat /= 0) CALL abort_program(routinen, "The XYZ file could not be opened")
397 IF (info) WRITE (unit=error_unit, fmt="(A)") "#", "# Reading XYZ file: "//trim(xyz_file_name)
398 READ (unit=xyz_input_unit, fmt="(A)", iostat=istat) arg
399 IF (istat /= 0) THEN
400 CALL abort_program(routinen, "Reading line 1 of the XYZ file (number of atoms) failed")
401 END IF
402 IF (arg(1:1) == "#") THEN
403 READ (unit=arg, fmt=*) string, natom_xyz
404 ELSE
405 READ (unit=arg, fmt=*) natom_xyz
406 END IF
407 IF (istat /= 0) THEN
408 CALL abort_program(routinen, "Reading line 1 of the XYZ file (number of atoms) failed")
409 END IF
410 IF (ALLOCATED(atomic_label)) DEALLOCATE (atomic_label)
411 ALLOCATE (atomic_label(natom_xyz), stat=istat)
412 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the vector atomic_label failed")
413 atomic_label(:) = ""
414 IF (ekin) THEN
415 IF (ALLOCATED(atomic_mass)) DEALLOCATE (atomic_mass)
416 ALLOCATE (atomic_mass(natom_xyz), stat=istat)
417 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the vector atomic_mass failed")
418 END IF
419 READ (unit=xyz_input_unit, fmt="(A)", iostat=istat) remark_xyz
420 IF (istat /= 0) CALL abort_program(routinen, "Reading line 2 of the XYZ file (remark line) failed")
421 iatom = 0
422 DO
423 READ (unit=xyz_input_unit, fmt=*) arg
424 IF (arg(1:1) == "#") THEN
425 cycle
426 ELSE
427 backspace(unit=xyz_input_unit)
428 END IF
429 iatom = iatom + 1
430 READ (unit=xyz_input_unit, fmt=*, iostat=istat) atomic_label(iatom), rdum(1:3)
431 IF (istat /= 0) THEN
432 message = ""
433 WRITE (unit=message, fmt="(A,I0,A,I0,A)") &
434 "Reading line ", iatom + 2, " of the XYZ file (atom ", iatom, ") failed"
435 CALL abort_program(routinen, trim(message))
436 END IF
437 CALL uppercase(atomic_label(iatom))
438 IF (len_trim(atomic_label(iatom)) > 1) CALL lowercase(atomic_label(iatom) (2:2))
439 atomic_label(iatom) = trim(adjustl(atomic_label(iatom)))
440 IF (ekin) atomic_mass(iatom) = get_atomic_mass(atomic_label(iatom))
441 IF (iatom == natom_xyz) EXIT
442 END DO
443 IF (info) THEN
444 WRITE (unit=error_unit, fmt="(A,I0)") &
445 "# Number of atoms : ", natom_xyz, &
446 "# Remark : "//trim(adjustl(remark_xyz))
447 END IF
448 END IF
449 END IF
450
451 ! Read cell information file if requested
452 IF (have_cell_file) THEN
453 string = ""
454 ! Check if a new cell file name was specified, if yes then close the old unit
455 INQUIRE (unit=cell_input_unit, name=string)
456 IF (trim(string) /= trim(cell_file_name)) CLOSE (unit=cell_input_unit)
457 INQUIRE (unit=cell_input_unit, opened=opened)
458 ! Read new cell file if needed and update reference information
459 IF (.NOT. opened) THEN
460 OPEN (unit=cell_input_unit, &
461 file=cell_file_name, &
462 status="OLD", &
463 access="SEQUENTIAL", &
464 form="FORMATTED", &
465 position="REWIND", &
466 action="READ", &
467 iostat=istat)
468 IF (istat /= 0) CALL abort_program(routinen, "The cell file could not be opened")
469 IF (info) WRITE (unit=error_unit, fmt="(A)") "#", "# Reading cell file: "//trim(cell_file_name)
470 END IF
471 END IF
472
473 ! Increment DCD file counter
474 ndcd_file = ndcd_file + 1
475
476 ! Open input file in DCD format
477 OPEN (unit=input_unit, &
478 file=dcd_file_name, &
479 status="OLD", &
480 form="UNFORMATTED", &
481 position="REWIND", &
482 action="READ", &
483 iostat=istat)
484 IF (istat /= 0) CALL abort_program(routinen, "The DCD file could not be opened")
485 IF (info) THEN
486 WRITE (unit=error_unit, fmt="(A,/,(A,I0))") &
487 "#", &
488 "# DCD file number : ", ndcd_file, &
489 "# Reading DCD file: "//trim(dcd_file_name)
490 END IF
491
492 ! Read 1st record of DCD file
493 READ (unit=input_unit) id_dcd, idum(1), istep_dcd, istride_dcd, idum(2:7), dt_dcd, have_unit_cell, idum(8:16)
494 IF (info) THEN
495 WRITE (unit=error_unit, fmt="(A,2(/,A,I0),/,A,F9.3,/,A,I0)") &
496 "# DCD id string : "//id_dcd, &
497 "# Step : ", istep_dcd, &
498 "# Print frequency : ", istride_dcd, &
499 "# Time step [fs] : ", dt_dcd, &
500 "# Unit cell : ", have_unit_cell
501 END IF
502
503 IF (ekin .AND. trim(adjustl(id_dcd)) /= "VEL") THEN
504 CALL abort_program(routinen, "ekin flag requires a DCD file with VELocities")
505 END IF
506 IF (apply_pbc .AND. (have_unit_cell /= 1)) THEN
507 CALL abort_program(routinen, "pbc flags require that unit cell information is available")
508 END IF
509
510 IF ((trim(adjustl(id_dcd)) == "FRC") .OR. (trim(adjustl(id_dcd)) == "VEL")) THEN
511 unit_string = "[a.u.]"
512 IF (apply_pbc) CALL abort_program(routinen, "pbc flags require a DCD file with COoRDinates")
513 ELSE IF (trim(adjustl(id_dcd)) == "CORD") THEN
514 unit_string = "[Angstrom]"
515 ELSE
516 CALL abort_program(routinen, "Unknown DCD id found (use -debug or -info flag for details)")
517 END IF
518
519 ! Read 2nd record of DCD file
520 READ (unit=input_unit) nremark, remark_dcd(1), remark_dcd(2)
521 IF (info) THEN
522 DO i = 1, nremark
523 WRITE (unit=error_unit, fmt="(A,I1,A)") &
524 "# Remark ", i, " : "//trim(remark_dcd(i))
525 END DO
526 END IF
527
528 ! Read 3rd record of DCD file
529 READ (unit=input_unit) natom_dcd
530 IF (info) THEN
531 WRITE (unit=error_unit, fmt="(A,I0)") &
532 "# Number of atoms : ", natom_dcd
533 END IF
534
535 IF (have_atomic_labels) THEN
536 IF (natom_dcd /= natom_xyz) CALL abort_program(routinen, "Number of atoms in XYZ and DCD file differ")
537 ELSE
538 IF (.NOT. ALLOCATED(atomic_label)) THEN
539 ALLOCATE (atomic_label(natom_dcd), stat=istat)
540 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the vector atomic_label failed")
541 atomic_label(:) = ""
542 END IF
543 END IF
544
545 ! Overwrite MD time step, if requested
546 IF (md_time_step > 0.0_dp) THEN
547 dt_dcd = real(md_time_step, kind=sp)
548 dt = md_time_step
549 ELSE
550 dt = real(dt_dcd, kind=dp)
551 END IF
552
553 ! Dump output in DCD format => dcd2dcd functionality
554 IF (output_format_dcd) THEN
555 IF (vel2cord) id_dcd = "CORD"
556 ! Note, dt_dcd is REAL*4 and the rest is INTEGER*4
557 WRITE (unit=output_unit) id_dcd, idum(1), istep_dcd, istride_dcd, idum(2:7), dt_dcd, have_unit_cell, idum(8:16)
558 WRITE (unit=output_unit) nremark, remark_dcd(1), remark_dcd(2)
559 WRITE (unit=output_unit) natom_dcd
560 END IF
561
562 frame_loop: DO
563
564 nframe = nframe + 1
565
566 IF (nframe < first_frame) THEN
567 dump_frame = .false.
568 ELSE
569 IF (modulo(nframe - first_frame, stride) == 0) THEN
570 dump_frame = .true.
571 ELSE
572 dump_frame = .false.
573 END IF
574 END IF
575
576 ! Read unit cell information, if available
577 IF (have_unit_cell == 1) THEN
578 READ (unit=input_unit, iostat=istat) a, gamma, b, beta, alpha, c
579 IF (istat < 0) EXIT frame_loop
580 END IF
581 IF (have_cell_file) THEN
582 ! Read cell information from cell file if available
583 DO
584 READ (unit=cell_input_unit, fmt=*, iostat=istat) string
585 IF (istat /= 0) THEN
586 message = ""
587 WRITE (unit=message, fmt="(A)") &
588 "Reading frame (step) from cell file "//trim(cell_file_name)//" failed"
589 CALL abort_program(routinen, trim(message))
590 END IF
591 IF (string(1:1) == "#") THEN
592 cycle
593 ELSE
594 backspace(unit=cell_input_unit)
595 EXIT
596 END IF
597 END DO
598 READ (unit=cell_input_unit, fmt=*, iostat=istat) iframe, step_time, hmat(1:3, 1:3), cell_volume
599 IF (istat /= 0) THEN
600 WRITE (unit=error_unit, fmt="(/,T2,A,I0)") &
601 "IOSTAT = ", istat, &
602 "Line read: """//trim(string)//""""
603 message = ""
604 WRITE (unit=message, fmt="(A)") &
605 "Invalid cell information read from cell file "//trim(cell_file_name)
606 CALL abort_program(routinen, trim(message))
607 END IF
608 ! Store step time of the first (selected) MD step
609 IF (nframe == first_frame) first_step_time = step_time
610 ! Perform checks only for the selected frames
611 IF ((nframe >= first_frame) .AND. (nframe <= last_frame)) THEN
612 ! Save cell information from DCD header, if available
613 IF (have_unit_cell == 1) THEN
614 a_dcd = a
615 b_dcd = b
616 c_dcd = c
617 alpha_dcd = alpha
618 beta_dcd = beta
619 gamma_dcd = gamma
620 END IF
621 ! Overwrite cell information from DCD header
622 a = norm2(hmat(1:3, 1))
623 b = norm2(hmat(1:3, 2))
624 c = norm2(hmat(1:3, 3))
625 ! Lattice vectors from DCD headers of VEL files are atomic units, but the CP2K cell file is in Angstrom
626 IF (unit_string == "[a.u.]") THEN
627 a = a/angstrom
628 b = b/angstrom
629 c = c/angstrom
630 END IF
631 alpha = angle(hmat(1:3, 2), hmat(1:3, 3))*degree
632 beta = angle(hmat(1:3, 3), hmat(1:3, 1))*degree
633 gamma = angle(hmat(1:3, 1), hmat(1:3, 2))*degree
634 IF (have_unit_cell == 1) THEN
635 ! Check consistency of DCD and cell file information
636 IF (.NOT. ignore_warnings) THEN
637 IF (md_time_step > 0.0_dp) THEN
638 string = "Time step (requested) ="
639 ELSE
640 string = "Time step (DCD header) ="
641 END IF
642 ndt = (step_time - first_step_time)/dt
643 IF (abs(ndt - anint(ndt)) > 0.01_dp) THEN
644 WRITE (unit=error_unit, fmt="(/,T2,A,I8,/,(T2,A,F15.6,A))") &
645 "Step number (CELL file) = ", iframe, &
646 "Step time (CELL file) = ", step_time, " fs", &
647 "First step (CELL file) = ", first_step_time, " fs", &
648 "Time since first step = ", step_time - first_step_time, " fs", &
649 "Steps since first step = ", ndt, "", &
650 trim(string)//" ", dt, " fs"
651 WRITE (unit=error_unit, fmt="(/,T2,A)") &
652 "*** WARNING: MD step time in cell file is not a multiple of the MD time step in the DCD file header ***"
653 END IF
654 END IF
655 eps_geo = 1.0e-7_dp
656 IF (abs(a - a_dcd) > eps_geo) THEN
657 WRITE (unit=error_unit, fmt="(/,(T2,A,F14.6))") &
658 "a (CELL file) = ", a, &
659 "a (DCD header) = ", a_dcd
660 CALL abort_program(routinen, "Cell file and DCD file information for lattice constant ""a"" differ")
661 END IF
662 IF (abs(b - b_dcd) > eps_geo) THEN
663 WRITE (unit=error_unit, fmt="(/,(T2,A,F14.6))") &
664 "b (CELL file) = ", b, &
665 "b (DCD header) = ", b_dcd
666 CALL abort_program(routinen, "Cell file and DCD file information for lattice constant ""b"" differ")
667 END IF
668 IF (abs(c - c_dcd) > eps_geo) THEN
669 WRITE (unit=error_unit, fmt="(/,(T2,A,F14.6))") &
670 "c (CELL file) = ", c, &
671 "c (DCD header) = ", c_dcd
672 CALL abort_program(routinen, "Cell file and DCD file information for lattice constant ""c"" differ")
673 END IF
674 eps_angle = 1.0e-4_dp
675 IF (abs(alpha - alpha_dcd) > eps_angle) THEN
676 WRITE (unit=error_unit, fmt="(/,(T2,A,F14.6))") &
677 "alpha (CELL file) = ", alpha, &
678 "alpha (DCD header) = ", alpha_dcd
679 CALL abort_program(routinen, "Cell file and DCD file information for cell angle ""alpha"" differ")
680 END IF
681 IF (abs(beta - beta_dcd) > eps_angle) THEN
682 WRITE (unit=error_unit, fmt="(/,(T2,A,F14.6))") &
683 "beta (CELL file) = ", beta, &
684 "beta (DCD header) = ", beta_dcd
685 CALL abort_program(routinen, "Cell file and DCD file information for cell angle ""beta"" differ")
686 END IF
687 IF (abs(gamma - gamma_dcd) > eps_angle) THEN
688 WRITE (unit=error_unit, fmt="(/,(T2,A,F14.6))") &
689 "gamma (CELL file) = ", gamma, &
690 "gamma (DCD header) = ", gamma_dcd
691 CALL abort_program(routinen, "Cell file and DCD file information for cell angle ""beta"" differ")
692 END IF
693 END IF
694 END IF
695 END IF
696
697 IF ((info .OR. trace_atoms) .AND. dump_frame) THEN
698 WRITE (unit=error_unit, fmt="(A,/,A,I0)") &
699 "#", "# Frame number : ", nframe
700 END IF
701
702 ! Print unit cell information, if available
703 IF (info .AND. dump_frame) THEN
704 IF (have_unit_cell == 1) THEN
705 WRITE (unit=error_unit, fmt="(A,/,(A,T19,A,F12.6))") &
706 "#", &
707 "# a "//trim(unit_string), ": ", a, &
708 "# b "//trim(unit_string), ": ", b, &
709 "# c "//trim(unit_string), ": ", c
710 WRITE (unit=error_unit, fmt="(A,F12.6)") &
711 "# alpha [degree] : ", alpha, &
712 "# beta [degree] : ", beta, &
713 "# gamma [degree] : ", gamma
714 END IF
715 END IF
716
717 ! Allocate the array for the current atomic positions if needed
718 IF (.NOT. ALLOCATED(r)) THEN
719 ALLOCATE (r(natom_dcd, 3), stat=istat)
720 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array r failed")
721 END IF
722
723 ! Read in the atomic positions of the current frame
724 READ (unit=input_unit, iostat=istat) r(1:natom_dcd, 1)
725 IF (istat < 0) EXIT frame_loop
726 READ (unit=input_unit) r(1:natom_dcd, 2)
727 READ (unit=input_unit) r(1:natom_dcd, 3)
728
729 ! Store the atomic positions of the first frame in the current DCD input file for
730 ! the output of the atomic displacements
731 IF ((nframe == 1) .AND. print_atomic_displacements) THEN
732 IF (.NOT. ALLOCATED(r0)) THEN
733 ALLOCATE (r0(natom_dcd, 3), stat=istat)
734 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array r0 failed")
735 END IF
736 r0(:, :) = r(:, :)
737 IF (.NOT. ALLOCATED(atomic_displacement)) THEN
738 ALLOCATE (atomic_displacement(natom_dcd), stat=istat)
739 IF (istat /= 0) THEN
740 CALL abort_program(routinen, "Allocation of the vector atomic_displacement failed")
741 END IF
742 END IF
743 atomic_displacement(:) = 0.0_dp
744 END IF
745
746 IF (ekin .AND. trim(adjustl(id_dcd)) == "VEL") THEN
747
748 ! Dump the kinetic energy of each as an "atomic temperature"
749 IF (dump_frame) THEN
750
751 IF (.NOT. ALLOCATED(atomic_temperature)) THEN
752 ALLOCATE (atomic_temperature(natom_dcd), stat=istat)
753 IF (istat /= 0) THEN
754 CALL abort_program(routinen, "Allocation of the vector atomic_temperature failed")
755 END IF
756 atomic_temperature(:) = 0.0_dp
757 END IF
758
759 IF (info) THEN
760 WRITE (unit=error_unit, fmt="(A)") &
761 "#", &
762 "# Temperature [K]"
763 END IF
764
765 tavg_frame = 0.0_dp
766 DO iatom = 1, natom_dcd
767 atomic_temperature(iatom) = atomic_mass(iatom)*(r(iatom, 1)*r(iatom, 1) + &
768 r(iatom, 2)*r(iatom, 2) + &
769 r(iatom, 3)*r(iatom, 3))*kelvin/3.0_dp
770 tavg_frame = tavg_frame + atomic_temperature(iatom)
771 END DO
772 tavg_frame = tavg_frame/real(natom_dcd, kind=dp)
773
774 IF (output_format_dcd) THEN
775 IF (have_unit_cell == 1) THEN
776 WRITE (unit=output_unit) a, gamma, b, beta, alpha, c
777 END IF
778 ! This is a hack for VMD: dump the atomic temperatures as the x-coordinates
779 ! of a DCD VEL file
780 WRITE (unit=output_unit) real(atomic_temperature(:), kind=sp)
781 ! The y and z coordinates are filled with zeros
782 atomic_temperature(:) = 0.0_dp
783 WRITE (unit=output_unit) real(atomic_temperature(:), kind=sp)
784 WRITE (unit=output_unit) real(atomic_temperature(:), kind=sp)
785 ELSE
786 DO iatom = 1, natom_dcd
787 WRITE (unit=output_unit, fmt="(A5,5X,F25.3)") adjustl(atomic_label(iatom)), atomic_temperature
788 END DO
789 END IF
790
791 IF (info) THEN
792 WRITE (unit=error_unit, fmt="(A,F12.3)") &
793 "# T [K] this frame: ", tavg_frame
794 END IF
795
796 tavg = tavg + tavg_frame
797
798 END IF ! dump_frame
799
800 ELSE
801
802 IF (dump_frame) THEN
803
804 ! Apply periodic boundary conditions before dumping the coordinates if requested
805 IF (apply_pbc) THEN
806 IF (.NOT. ALLOCATED(r_pbc)) THEN
807 ALLOCATE (r_pbc(natom_dcd, 3), stat=istat)
808 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array r_pbc failed")
809 r_pbc(:, :) = 0.0_dp
810 END IF
811 IF (.NOT. ALLOCATED(s)) THEN
812 ALLOCATE (s(natom_dcd, 3), stat=istat)
813 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array s failed")
814 s(:, :) = 0.0_dp
815 END IF
816 IF (.NOT. ALLOCATED(s_pbc)) THEN
817 ALLOCATE (s_pbc(natom_dcd, 3), stat=istat)
818 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array s_pbc failed")
819 s_pbc(:, :) = 0.0_dp
820 END IF
821 CALL pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
822 CALL write_out_of_box_atoms(atomic_label, r, s, eps_out_of_box, h)
823 ! Overwrite input coordinate with the PBCed coordinates for printing
824 r(:, :) = real(r_pbc(:, :), kind=sp)
825 END IF ! apply_pbc
826
827 ! Calculate the atomic displacements with respect to the first frame,
828 ! i.e. set of atomic positions
829 IF (print_atomic_displacements) THEN
830 DO iatom = 1, natom_dcd
831 atomic_displacement(iatom) = sqrt((r(iatom, 1) - r0(iatom, 1))**2 + &
832 (r(iatom, 2) - r0(iatom, 2))**2 + &
833 (r(iatom, 3) - r0(iatom, 3))**2)
834 END DO
835 END IF
836
837 ! Dump the atomic coordinates in DCD or XMOL/MOLDEN format
838 IF (output_format_dcd) THEN
839 IF (have_unit_cell == 1) THEN
840 WRITE (unit=output_unit) a, gamma, b, beta, alpha, c
841 END IF
842 IF (print_atomic_displacements) THEN
843 ! This is a hack for VMD: dump the atomic displacements as the x-coordinates
844 ! of a DCD CORD file
845 WRITE (unit=output_unit) real(atomic_displacement(:), kind=sp)
846 ! The y and z coordinates are filled with zeros
847 atomic_displacement(:) = 0.0_dp
848 WRITE (unit=output_unit) real(atomic_displacement(:), kind=sp)
849 WRITE (unit=output_unit) real(atomic_displacement(:), kind=sp)
850 ELSE
851 ! Dump DCD file information
852 DO i = 1, 3
853 WRITE (unit=output_unit) r(:, i)
854 END DO
855 END IF
856 ELSE
857 IF (print_atomic_displacements) THEN
858 IF (info) THEN
859 WRITE (unit=error_unit, fmt="(A,2(/,A),T15,A)") &
860 "#", &
861 "# Displacements :", &
862 "# "//trim(unit_string), "|r|"
863 END IF
864 IF (have_atomic_labels) THEN
865 DO iatom = 1, natom_dcd
866 WRITE (unit=output_unit, fmt="(A5,5X,F25.6)") adjustl(atomic_label(iatom)), atomic_displacement(iatom)
867 END DO
868 ELSE
869 DO iatom = 1, natom_dcd
870 WRITE (unit=output_unit, fmt="(10X,F25.6)") atomic_displacement(iatom)
871 END DO
872 END IF
873 ELSE
874 IF (output_format_xmol) THEN
875 IF (have_cell_file) THEN
876 WRITE (unit=output_unit, fmt="(T2,I0,/,3(1X,F14.6),3(1X,F9.3))") &
877 natom_dcd, a, b, c, alpha, beta, gamma
878 ELSE
879 IF (have_unit_cell == 1) THEN
880 WRITE (unit=output_unit, fmt="(T2,I0,/,A,I8,A,F12.3,A,F20.10,3(A,F14.6),3(A,F8.3))") &
881 natom_dcd, " i = ", istep_dcd + nframe - 1, &
882 ", time = ", real((nframe - 1)**istride_dcd, kind=dp)*dt, &
883 ", E = ", energy, &
884 ", a = ", a, &
885 ", b = ", b, &
886 ", c = ", c, &
887 ", alpha = ", alpha, &
888 ", beta = ", beta, &
889 ", gamma = ", gamma
890 ELSE
891 WRITE (unit=output_unit, fmt="(T2,I0,/,A,I8,A,F12.3,A,F20.10)") &
892 natom_dcd, " i = ", istep_dcd + nframe - 1, &
893 ", time = ", real((nframe - 1)**istride_dcd, kind=dp)*dt, &
894 ", E = ", energy
895 END IF
896 END IF
897 DO iatom = 1, natom_dcd
898 WRITE (unit=output_unit, fmt=fmt_string) adjustl(atomic_label(iatom)), r(iatom, 1:3)
899 END DO
900 ELSE
901 IF (info) THEN
902 WRITE (unit=error_unit, fmt="(A,T20,A)") &
903 "# "//trim(unit_string), "x y z"
904 END IF
905 IF (print_scaled_coordinates) THEN
906 DO iatom = 1, natom_dcd
907 WRITE (unit=output_unit, fmt=fmt_string) adjustl(atomic_label(iatom)), s(iatom, 1:3)
908 END DO
909 ELSE IF (print_scaled_pbc_coordinates) THEN
910 DO iatom = 1, natom_dcd
911 WRITE (unit=output_unit, fmt=fmt_string) adjustl(atomic_label(iatom)), s_pbc(iatom, 1:3)
912 END DO
913 ELSE
914 DO iatom = 1, natom_dcd
915 WRITE (unit=output_unit, fmt=fmt_string) adjustl(atomic_label(iatom)), r(iatom, 1:3)
916 END DO
917 END IF
918 END IF
919 END IF
920 END IF ! output_format_dcd
921
922 END IF ! dump_frame
923
924 END IF
925
926 IF (dump_frame) nframe_read = nframe_read + 1
927
928 ! Exit loop and stop processing, if the last (requested) frame was encountered
929 IF (nframe >= last_frame) THEN
930 nframe = nframe - 1
931 CLOSE (unit=input_unit)
932 EXIT dcd_file_loop
933 END IF
934
935 END DO frame_loop
936
937 nframe = nframe - 1
938
939 CLOSE (unit=input_unit)
940
941 IF (iarg >= narg) EXIT dcd_file_loop
942
943 END DO dcd_file_loop
944
945 IF (info) THEN
946 WRITE (unit=error_unit, fmt="(A,/,A,I0)") &
947 "#", &
948 "# Frames processed: ", nframe_read
949 END IF
950
951 IF (ekin .AND. trim(adjustl(id_dcd)) == "VEL") THEN
952 IF (info) THEN
953 WRITE (unit=error_unit, fmt="(A,/,A,F12.3)") &
954 "#", &
955 "# T [K] all frames: ", tavg/real(nframe_read, kind=dp)
956 END IF
957 END IF
958
959 ! Print version information
960 IF (info) THEN
961 WRITE (unit=error_unit, fmt="(A)") &
962 "#", &
963 "# Normal termination of "//trim(version_info)
964 END IF
965
966 ! Close files
967 IF (have_atomic_labels) CLOSE (unit=xyz_input_unit)
968 IF (len_trim(out_file_name) > 0) CLOSE (unit=output_unit)
969
970 ! Cleanup
971 IF (ALLOCATED(atomic_label)) DEALLOCATE (atomic_label)
972 IF (ALLOCATED(atomic_displacement)) DEALLOCATE (atomic_displacement)
973 IF (ALLOCATED(atomic_mass)) DEALLOCATE (atomic_mass)
974 IF (ALLOCATED(atomic_temperature)) DEALLOCATE (atomic_temperature)
975 IF (ALLOCATED(r)) DEALLOCATE (r)
976 IF (ALLOCATED(r0)) DEALLOCATE (r0)
977 IF (ALLOCATED(r_pbc)) DEALLOCATE (r_pbc)
978 IF (ALLOCATED(s)) DEALLOCATE (s)
979 IF (ALLOCATED(s_pbc)) DEALLOCATE (s_pbc)
980
981CONTAINS
982
983! **************************************************************************************************
984!> \brief ...
985!> \param routine ...
986!> \param message ...
987! **************************************************************************************************
988 SUBROUTINE abort_program(routine, message)
989 ! Abort the program after printing an error message to stderr
990
991 CHARACTER(LEN=*), INTENT(IN) :: routine, message
992
993 CHARACTER(LEN=2*default_string_length) :: error_message
994
995 error_message = "*** ERROR in "//trim(routine)//": "//trim(message)//" ***"
996 WRITE (unit=default_error_unit, fmt="(/,A,/)") trim(error_message)
997 stop "*** ABNORMAL PROGRAM TERMINATION of dumpdcd v3.2 ***"
998
999 END SUBROUTINE abort_program
1000
1001! **************************************************************************************************
1002!> \brief Calculation of the angle between the vectors a and b.
1003!> The angle is returned in radians.
1004!> \param a ...
1005!> \param b ...
1006!> \return ...
1007! **************************************************************************************************
1008 PURE FUNCTION angle(a, b) RESULT(angle_ab)
1009
1010 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: a, b
1011 REAL(kind=dp) :: angle_ab
1012
1013 REAL(kind=dp), PARAMETER :: eps_geo = 1.0e-6_dp
1014
1015 REAL(kind=dp) :: length_of_a, length_of_b
1016 REAL(kind=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
1017
1018 length_of_a = norm2(a)
1019 length_of_b = norm2(b)
1020
1021 IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
1022 a_norm(:) = a(:)/length_of_a
1023 b_norm(:) = b(:)/length_of_b
1024 angle_ab = acos(min(max(dot_product(a_norm, b_norm), -1.0_dp), 1.0_dp))
1025 ELSE
1026 angle_ab = 0.0_dp
1027 END IF
1028
1029 END FUNCTION angle
1030
1031! **************************************************************************************************
1032!> \brief ...
1033!> \param a ...
1034!> \param b ...
1035!> \param c ...
1036!> \param alpha ...
1037!> \param beta ...
1038!> \param gamma ...
1039!> \param h ...
1040! **************************************************************************************************
1041 SUBROUTINE build_h_matrix(a, b, c, alpha, beta, gamma, h)
1042 ! Calculate the h matrix of a simulation cell given the cell edge lengths a, b,
1043 ! and c in Angstrom and the angles alpha (<b,c)), beta (<(a,c)), and gamma (<(a,b))
1044 ! in degree.
1045
1046 REAL(kind=dp), INTENT(IN) :: a, b, c, alpha, beta, gamma
1047 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: h
1048
1049 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_h_matrix'
1050
1051 REAL(kind=dp) :: cosa, cosb, cosg, sing
1052
1053 cosa = cos(alpha/degree)
1054 IF (abs(cosa) < epsilon(0.0_dp)) cosa = 0.0_dp
1055
1056 cosb = cos(beta/degree)
1057 IF (abs(cosb) < epsilon(0.0_dp)) cosb = 0.0_dp
1058
1059 cosg = cos(gamma/degree)
1060 IF (abs(cosg) < epsilon(0.0_dp)) cosg = 0.0_dp
1061
1062 sing = sin(gamma/degree)
1063 IF (abs(sing) < epsilon(0.0_dp)) sing = 0.0_dp
1064
1065 h(1, 1) = 1.0_dp
1066 h(2, 1) = 0.0_dp
1067 h(3, 1) = 0.0_dp
1068
1069 h(1, 2) = cosg
1070 h(2, 2) = sing
1071 h(3, 2) = 0.0_dp
1072
1073 h(1, 3) = cosb
1074 h(2, 3) = (cosa - cosg*cosb)/sing
1075 IF ((1.0_dp - h(1, 3)**2 - h(2, 3)**2) < 0.0_dp) THEN
1076 CALL abort_program(routinen, "Build of the h matrix failed, check cell information")
1077 END IF
1078 h(3, 3) = sqrt(1.0_dp - h(1, 3)**2 - h(2, 3)**2)
1079
1080 h(:, 1) = a*h(:, 1)
1081 h(:, 2) = b*h(:, 2)
1082 h(:, 3) = c*h(:, 3)
1083
1084 END SUBROUTINE build_h_matrix
1085
1086! **************************************************************************************************
1087!> \brief ...
1088!> \param a ...
1089!> \return ...
1090! **************************************************************************************************
1091 FUNCTION det_3x3(a) RESULT(det_a)
1092 ! Returns the determinante of the 3x3 matrix a.
1093
1094 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: a
1095 REAL(kind=dp) :: det_a
1096
1097 det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
1098 a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
1099 a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
1100
1101 END FUNCTION det_3x3
1102
1103! **************************************************************************************************
1104!> \brief ...
1105!> \param element_symbol ...
1106!> \return ...
1107! **************************************************************************************************
1108 FUNCTION get_atomic_mass(element_symbol) RESULT(amass)
1109 ! Get the atomic mass amass for an element
1110
1111 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
1112 REAL(kind=dp) :: amass
1113
1114 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_atomic_mass'
1115
1116 SELECT CASE (trim(element_symbol))
1117 CASE ("O ")
1118 amass = 15.9994_dp
1119 CASE ("U ")
1120 amass = 238.02891_dp
1121 CASE DEFAULT
1122 CALL abort_program(routinen, "Unknown element symbol found")
1123 END SELECT
1124
1125 amass = amass*massunit
1126
1127 END FUNCTION get_atomic_mass
1128
1129! **************************************************************************************************
1130!> \brief ...
1131!> \param h ...
1132!> \param hinv ...
1133!> \param deth ...
1134! **************************************************************************************************
1135 SUBROUTINE invert_matrix_3x3(h, hinv, deth)
1136 ! Calculate the inverse hinv and the determinant deth of the 3x3 matrix h.
1137
1138 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: h
1139 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: hinv
1140 REAL(kind=dp), INTENT(OUT) :: deth
1141
1142 CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_matrix_3x3'
1143
1144 deth = det_3x3(h)
1145
1146 ! Numerics
1147 deth = abs(deth)
1148 IF (deth < 1.0e-10_dp) THEN
1149 CALL abort_program(routinen, "Invalid h matrix for cell found; det(h) < 1.0E-10")
1150 END IF
1151
1152 hinv(1, 1) = (h(2, 2)*h(3, 3) - h(3, 2)*h(2, 3))/deth
1153 hinv(2, 1) = (h(2, 3)*h(3, 1) - h(3, 3)*h(2, 1))/deth
1154 hinv(3, 1) = (h(2, 1)*h(3, 2) - h(3, 1)*h(2, 2))/deth
1155
1156 hinv(1, 2) = (h(1, 3)*h(3, 2) - h(3, 3)*h(1, 2))/deth
1157 hinv(2, 2) = (h(1, 1)*h(3, 3) - h(3, 1)*h(1, 3))/deth
1158 hinv(3, 2) = (h(1, 2)*h(3, 1) - h(3, 2)*h(1, 1))/deth
1159
1160 hinv(1, 3) = (h(1, 2)*h(2, 3) - h(2, 2)*h(1, 3))/deth
1161 hinv(2, 3) = (h(1, 3)*h(2, 1) - h(2, 3)*h(1, 1))/deth
1162 hinv(3, 3) = (h(1, 1)*h(2, 2) - h(2, 1)*h(1, 2))/deth
1163
1164 END SUBROUTINE invert_matrix_3x3
1165
1166! **************************************************************************************************
1167!> \brief ...
1168!> \param string ...
1169! **************************************************************************************************
1170 SUBROUTINE lowercase(string)
1171 ! Convert all letters in a string to lowercase
1172 CHARACTER(LEN=*), INTENT(INOUT) :: string
1173
1174 INTEGER :: i, iascii
1175
1176 DO i = 1, len_trim(string)
1177 iascii = ichar(string(i:i))
1178 IF ((iascii >= 65) .AND. (iascii <= 90)) THEN
1179 string(i:i) = char(iascii + 32)
1180 END IF
1181 END DO
1182
1183 END SUBROUTINE lowercase
1184
1185! **************************************************************************************************
1186!> \brief ...
1187!> \param r ...
1188!> \param r_pbc ...
1189!> \param s ...
1190!> \param s_pbc ...
1191!> \param a ...
1192!> \param b ...
1193!> \param c ...
1194!> \param alpha ...
1195!> \param beta ...
1196!> \param gamma ...
1197!> \param debug ...
1198!> \param info ...
1199!> \param pbc0 ...
1200!> \param h ...
1201!> \param hinv ...
1202! **************************************************************************************************
1203 SUBROUTINE pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
1204 ! Apply the periodic boundary conditions (PBC) to the Cartesian coordinate array
1205 ! r given the cell edge lengths a, b, and c in Angstrom and the angles alpha (<b,c)),
1206 ! beta (<(a,c)), and gamma (<(a,b)) in degree.
1207 ! On output r_pbc is updated with the "PBCed" input coordinates, s with the scaled
1208 ! input coordinates, and s_pbc with scaled "PBCed" coordinates.
1209 ! If pbc0 is true then fold to the range [-l/2,+l/2[ (origin at box centre) else fold
1210 ! to the range [0,l[ (origin at lower left box corner).
1211
1212 REAL(kind=sp), DIMENSION(:, :), INTENT(IN) :: r
1213 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: r_pbc, s, s_pbc
1214 REAL(kind=dp), INTENT(IN) :: a, b, c, alpha, beta, gamma
1215 LOGICAL, INTENT(IN) :: debug, info, pbc0
1216 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: h, hinv
1217
1218 CHARACTER(LEN=*), PARAMETER :: routinen = 'pbc'
1219
1220 INTEGER :: i, natom
1221 LOGICAL :: orthorhombic
1222 REAL(kind=dp) :: deth
1223
1224 natom = SIZE(r, 1)
1225 IF (SIZE(r, 2) /= 3) CALL abort_program(routinen, "Array dimension for r must be 3")
1226
1227 ! Build h matrix
1228
1229 h(:, :) = 0.0_dp
1230
1231 IF ((abs(alpha - 90.0_dp) < epsilon(0.0_dp)) .AND. &
1232 (abs(beta - 90.0_dp) < epsilon(0.0_dp)) .AND. &
1233 (abs(gamma - 90.0_dp) < epsilon(0.0_dp))) THEN
1234 orthorhombic = .true.
1235 h(1, 1) = a
1236 h(2, 2) = b
1237 h(3, 3) = c
1238 ELSE
1239 orthorhombic = .false.
1240 CALL build_h_matrix(a, b, c, alpha, beta, gamma, h)
1241 END IF
1242
1243 ! Build inverse of h
1244 hinv(:, :) = 0.0_dp
1245 CALL invert_matrix_3x3(h, hinv, deth)
1246
1247 IF (info) THEN
1248 WRITE (unit=error_unit, fmt="(A)") "#"
1249 IF (orthorhombic) THEN
1250 WRITE (unit=error_unit, fmt="(A)") "# Cell symmetry : orthorhombic"
1251 ELSE
1252 WRITE (unit=error_unit, fmt="(A)") "# Cell symmetry : non-orthorhombic"
1253 END IF
1254 IF (debug) THEN
1255 WRITE (unit=error_unit, fmt="(A)") "#"
1256 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# / ", h(1, :), " \"
1257 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# h = | ", h(2, :), " |"
1258 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# \ ", h(3, :), " /"
1259 WRITE (unit=error_unit, fmt="(A)") "#"
1260 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# / ", hinv(1, :), " \"
1261 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# Inv(h) = | ", hinv(2, :), " |"
1262 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# \ ", hinv(3, :), " /"
1263 WRITE (unit=error_unit, fmt="(A)") "#"
1264 WRITE (unit=error_unit, fmt="(A,F0.6)") "# det(h) = ", deth
1265 END IF
1266 END IF
1267
1268 ! Calculate scaled coordinates and wrap back all atoms, i.e. apply the PBC
1269 IF (orthorhombic) THEN
1270 ! Try to save some flops in the case of an orthorhombic box
1271 DO i = 1, 3
1272 s(:, i) = r(:, i)*hinv(i, i)
1273 END DO
1274 ELSE
1275 s(:, :) = matmul(r(:, :), transpose(hinv(:, :)))
1276 END IF
1277
1278 IF (pbc0) THEN
1279 s_pbc(:, :) = s(:, :) - anint(s(:, :))
1280 ELSE
1281 s_pbc(:, :) = s(:, :) - floor(s(:, :))
1282 END IF
1283
1284 IF (orthorhombic) THEN
1285 DO i = 1, 3
1286 r_pbc(:, i) = s_pbc(:, i)*h(i, i)
1287 END DO
1288 ELSE
1289 r_pbc(:, :) = matmul(s_pbc(:, :), transpose(h(:, :)))
1290 END IF
1291
1292 END SUBROUTINE pbc
1293
1294! **************************************************************************************************
1295!> \brief ...
1296! **************************************************************************************************
1297 SUBROUTINE print_help()
1298
1299 ! Print the program flags for help
1300
1301 WRITE (unit=*, fmt="(T2,A)") &
1302 "", &
1303 "Program flags for "//trim(version_info)//":", &
1304 "", &
1305 " -cell_file, -cell : Input file with cell information in CP2K format (.cell)", &
1306 " -debug, -d : Print debug information", &
1307 " -ekin : Dump just the ""temperature"" of each atom", &
1308 " -eformat : Print coordinates in scientific format", &
1309 " -eo : Write standard output and standard error to the same logical unit", &
1310 " -first_frame, -ff <int> : Number of the first frame which is dumped", &
1311 " -help, -h : Print this information", &
1312 " -ignore_warnings : Do not print warning messages, e.g. about inconsistencies", &
1313 " -info, -i : Print additional information for each frame (see also -debug flag)", &
1314 " -last_frame, -lf <int> : Number of the last frame which is dumped", &
1315 " -md_time_step <real> : Use this MD time step instead of the one from the DCD header", &
1316 " -output, -o <file_name> : Name of the output file (default is stdout)", &
1317 " -output_format, -of <DCD|XMOL> : Output format for dump", &
1318 " -pbc : Apply the periodic boundary conditions (PBC) to each frame before it is dumped", &
1319 " (origin at lower left)", &
1320 " -pbc0 : Apply the periodic boundary conditions (PBC) to each frame before it is dumped", &
1321 " (origin at box centre)", &
1322 " -stride <int> : Stride for frame dump (allows to skip frames, e.g. by dumping each 10th frame)", &
1323 " -trace_atoms <real> : Print atoms which left the simulation box given a threshold value in scaled units", &
1324 " -vel2cord, -v2c : Dump a VELocity DCD file as COoRDinate DCD file (hack which allows a digest by VMD)", &
1325 " -xyz_file, -xyz <file_name> : Name of a reference XYZ file in XMOL format that provides the atomic labels", &
1326 ""
1327
1328 WRITE (unit=*, fmt="(T2,A)") &
1329 "Usage examples:", &
1330 "", &
1331 " dumpdcd <optional flags> <DCD file(s)>", &
1332 "", &
1333 "Specific usage examples:", &
1334 "", &
1335 " dumpdcd project-pos-1.dcd (without atomic labels from XYZ file)", &
1336 " dumpdcd -xyz project.xyz project-pos-1.dcd (single DCD file)", &
1337 " dumpdcd -xyz project.xyz project-pos-1.dcd project-pos-2.dcd ... (multiple DCD files are dumped consecutively)", &
1338 " dumpdcd -xyz project.xyz -cell_file project-1.cell -of xmol project-pos-1.dcd (check and print cell parameters in ", &
1339 " XMOL comment line, e.g. for TRAVIS, instead of the default REFTRAJ information)", &
1340 " dumpdcd -info -xyz project.xyz project-pos-1.dcd project-pos-2.dcd (print additional information)", &
1341 " dumpdcd -debug -xyz project.xyz project-pos-1.dcd project-pos-2.dcd (print debug information)", &
1342 " dumpdcd -ekin -d -xyz project.xyz project-vel-1.dcd (print the ""temperature"" of each atom)", &
1343 " dumpdcd -ekin -xyz project.xyz project-vel-1.dcd (print just the temperature of each atom)", &
1344 " dumpdcd -first_frame 5 -last_frame 10 project-pos-1.dcd (just dump frame 5 to 10, ie. 6 frames in total)", &
1345 " dumpdcd -o outfile.xyz project-pos-1.dcd (write output to the file ""outfile.xyz"" instead of stdout)", &
1346 " dumpdcd -o test.xyz -output_format xmol -xyz ref.xyz -first 10 -last 10 test.dcd (dump 10th frame in XMOL format)", &
1347 " dumpdcd -of dcd -ff 10 -lf 20 test.dcd (dump the frames 10 to 20 in DCD format to the default output file output.dcd)", &
1348 " dumpdcd -o part.dcd -of dcd -ff 1 -lf 3 test.dcd (dump the frames 1 to 3 in DCD format to the output file part.dcd)", &
1349 " dumpdcd -o part.dcd -of dcd -first 10 -lf 100 -stride 10 test.dcd (dump the frames 10,..., 100 to the file part.dcd)", &
1350 " dumpdcd -output new.dcd -output_format dcd -pbc old.dcd (dump all frames applying PBC to the output file new.dcd)", &
1351 " dumpdcd -o new.dcd -of dcd -pbc -trace_atoms 0.02 old.dcd (all atoms more than 2% out of the box are listed)", &
1352 " dumpdcd -o new.dcd -e out_of_box.log -of dcd -pbc -trace_atoms 0.1 old.dcd (atoms >10% out of the box are listed)", &
1353 " dumpdcd -o new.dcd -of dcd -vel2cord old.dcd (dump old.dcd as new.dcd and change only the DCD id from VEL to CORD)", &
1354 " dumpdcd -o new.dcd -of dcd -frc2cord old.dcd (dump old.dcd as new.dcd and change only the DCD id from FRC to CORD)", &
1355 " dumpdcd -i -disp project-pos-1.dcd (dump the displacements of all atoms w.r.t. their positions in the first frame)", &
1356 " dumpdcd -i -of dcd -disp project-pos-1.dcd (dump the atomic displacements as x-coordinates of a DCD CORD file)", &
1357 " dumpdcd -i -of dcd -ekin -v2c -xyz project.xyz project-vel-1.dcd (dump the atomic temperatures as x-coordinates of a ", &
1358 " DCD CORD file -> hack for VMD)", &
1359 ""
1360
1361 WRITE (unit=*, fmt="(T2,A)") &
1362 "Notes:", &
1363 "", &
1364 " - For -ekin a XYZ file is required to obtain the atomic labels", &
1365 " - The -info and the -debug flags provide a more detailed output which is especially handy for tracing problems", &
1366 " - The output in DCD format is in binary format", &
1367 " - The input coordinates should be in Angstrom. Velocities and forces are expected to be in atomic units", &
1368 ""
1369
1370 END SUBROUTINE print_help
1371
1372! **************************************************************************************************
1373!> \brief ...
1374!> \param string ...
1375! **************************************************************************************************
1376 SUBROUTINE uppercase(string)
1377 ! Convert all letters in a string to uppercase
1378
1379 CHARACTER(LEN=*), INTENT(INOUT) :: string
1380
1381 INTEGER :: i, iascii
1382
1383 DO i = 1, len_trim(string)
1384 iascii = ichar(string(i:i))
1385 IF ((iascii >= 97) .AND. (iascii <= 122)) THEN
1386 string(i:i) = char(iascii - 32)
1387 END IF
1388 END DO
1389
1390 END SUBROUTINE uppercase
1391
1392! **************************************************************************************************
1393!> \brief ...
1394!> \param atomic_label ...
1395!> \param r ...
1396!> \param s ...
1397!> \param eps_out_of_box ...
1398!> \param h ...
1399! **************************************************************************************************
1400 SUBROUTINE write_out_of_box_atoms(atomic_label, r, s, eps_out_of_box, h)
1401 ! Print a list of all atoms which have left the simulation box
1402
1403 CHARACTER(LEN=5), DIMENSION(:), INTENT(IN) :: atomic_label
1404 REAL(kind=sp), DIMENSION(:, :), INTENT(IN) :: r
1405 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: s
1406 REAL(kind=dp), INTENT(IN) :: eps_out_of_box
1407 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: h
1408
1409 INTEGER :: i, iatom, natom, ncount
1410 REAL(kind=dp) :: rl, s_max, s_min, sl
1411 REAL(kind=dp), DIMENSION(3) :: dr, ds
1412
1413 ! Quick return, if no action is requested
1414 IF (eps_out_of_box <= 0.0_dp) RETURN
1415
1416 s_max = 1.0_dp + eps_out_of_box
1417 s_min = -eps_out_of_box
1418 natom = SIZE(s, 1)
1419 ncount = 0
1420 DO iatom = 1, natom
1421 IF (any(s(iatom, :) < s_min) .OR. &
1422 any(s(iatom, :) > s_max)) THEN
1423 ncount = ncount + 1
1424 IF (ncount == 1) THEN
1425 WRITE (unit=error_unit, fmt="(A)") &
1426 "#", &
1427 "# Atoms out of box:", &
1428 "# Atom index label x y z |dr| |ds|"
1429 END IF
1430 ds(:) = s(iatom, :)
1431 DO i = 1, 3
1432 IF (s(iatom, i) < 0.0_dp) ds(i) = 0.0_dp
1433 IF (s(iatom, i) >= 1.0_dp) ds(i) = 1.0_dp
1434 END DO
1435 ds(:) = s(iatom, :) - ds(:)
1436 sl = sqrt(ds(1)**2 + ds(2)**2 + ds(3)**2)
1437 dr(:) = matmul(h(:, :), ds(:))
1438 rl = sqrt(dr(1)**2 + dr(2)**2 + dr(3)**2)
1439 WRITE (unit=error_unit, fmt="(A,I10,1X,A5,5(1X,F14.6))") &
1440 "# ", iatom, adjustr(atomic_label(iatom)), r(iatom, :), rl, sl
1441 END IF
1442 END DO
1443 WRITE (unit=error_unit, fmt="(A,I0,A)") "# ", ncount, " atom(s) out of box"
1444
1445 END SUBROUTINE write_out_of_box_atoms
1446
1447END PROGRAM dumpdcd
dumpdcd
program dumpdcd
Definition dumpdcd.F:8
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15