(git:e7e05ae)
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 
8 PROGRAM 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
192 100 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
212 101 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
222 102 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
264 104 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
275 108 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 
980 CONTAINS
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 
1446 END PROGRAM dumpdcd
abort_program
subroutine abort_program(routine, message)
...
Definition: dumpdcd.F:988
write_out_of_box_atoms
subroutine write_out_of_box_atoms(atomic_label, r, s, eps_out_of_box, h)
...
Definition: dumpdcd.F:1400
invert_matrix_3x3
subroutine invert_matrix_3x3(h, hinv, deth)
...
Definition: dumpdcd.F:1135
build_h_matrix
subroutine build_h_matrix(a, b, c, alpha, beta, gamma, h)
...
Definition: dumpdcd.F:1041
det_3x3
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
print_help
subroutine print_help()
...
Definition: dumpdcd.F:1297
get_atomic_mass
real(kind=dp) function get_atomic_mass(element_symbol)
...
Definition: dumpdcd.F:1108
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
lowercase
subroutine lowercase(string)
...
Definition: dumpdcd.F:1170
angle
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: dumpdcd.F:1008
uppercase
subroutine uppercase(string)
...
Definition: dumpdcd.F:1376
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