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