(git:71c3ab0)
Loading...
Searching...
No Matches
xyz2dcd.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8PROGRAM xyz2dcd
9
10! Version: 1.0
11! Author: Matthias Krack (MK)
12! History: - Creation (30.03.2015,MK)
13!
14! Note: The input coordinates and the cell vectors should be in Angstrom.
15
16! Uncomment the following line if this module is available (e.g. with gfortran)
17! and comment the corresponding variable declarations below
18! USE ISO_FORTRAN_ENV, ONLY: error_unit,input_unit,output_unit
19
20 IMPLICIT NONE
21
22 ! Comment the following lines if the ISO_FORTRAN_ENV is used (see above)
23 INTEGER, PARAMETER :: default_error_unit = 0, &
24 default_input_unit = 5, &
25 default_output_unit = 6
26 INTEGER :: error_unit = default_error_unit, &
27 output_unit = default_output_unit
28 ! End Comment
29
30 ! Parameters
31 CHARACTER(LEN=*), PARAMETER :: routinen = "xyz2dcd", &
32 version_info = routinen//" v1.0 (30.06.2020, Matthias Krack)"
33
34 INTEGER, PARAMETER :: dp = selected_real_kind(14, 200), &
35 sp = selected_real_kind(6, 30)
36 INTEGER, PARAMETER :: default_string_length = 240, &
37 cell_file_unit = 10, &
38 dcd_file_unit = 11, &
39 xyz_file_unit = default_input_unit
40
41 REAL(kind=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
42! REAL(KIND=dp), PARAMETER :: angstrom = 0.52917720859_dp ! [a.u.] -> [Angstrom]
43 REAL(kind=dp), PARAMETER :: degree = 180.0_dp/pi ! [rad] -> [degree]
44
45 ! Variables
46 CHARACTER(LEN=default_string_length) :: arg, cell_file_name, dcd_file_name, message, remark1, remark2, remark_xyz, &
47 string, xyz_file_name
48 CHARACTER(LEN=5), DIMENSION(:), ALLOCATABLE :: atomic_label
49 CHARACTER(LEN=80), DIMENSION(2) :: remark_dcd
50 INTEGER :: first_frame, i, iarg, iatom, iskip, istat, istep, j, &
51 last_frame, narg, natom, nframe, nframe_read, nremark, stride
52 LOGICAL :: apply_pbc, debug, dump_frame, exists, have_cell_file, have_cell_info, &
53 info, pbc0, print_atomic_displacements, print_scaled_coordinates, &
54 print_scaled_pbc_coordinates, trace_atoms
55 REAL(kind=dp) :: alpha, beta, dt, eps_out_of_box, gamma, tstep
56 REAL(kind=dp), DIMENSION(3) :: a, abc, b, c
57 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: atomic_displacement
58 REAL(kind=dp), DIMENSION(3, 3) :: h, hinv
59 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: r, r_pbc, r0, s, s_pbc
60
61 apply_pbc = .false.
62 debug = .false.
63 dump_frame = .true.
64 have_cell_file = .false.
65 have_cell_info = .false.
66 info = .false.
67 trace_atoms = .false.
68 pbc0 = .false.
69 print_scaled_coordinates = .false.
70 print_atomic_displacements = .false.
71 print_scaled_pbc_coordinates = .false.
72 first_frame = 1
73 last_frame = 1000000 ! Hard limit of 1 Mio frames in total
74 stride = 1
75 nframe = 0
76 nframe_read = 0
77 nremark = 0
78 cell_file_name = ""
79 dcd_file_name = ""
80 xyz_file_name = ""
81 remark_dcd(:) = ""
82 remark_xyz = ""
83 eps_out_of_box = -huge(0.0_dp)
84 h(:, :) = 0.0_dp
85 hinv(:, :) = 0.0_dp
86
87 ! Scan argument list and digest it
88
89 narg = command_argument_count()
90
91 IF (narg == 0) THEN
92 CALL print_help()
93 CALL abort_program(routinen, "No input file(s) specified")
94 END IF
95
96 iarg = 0
97 arg_loop: DO
98 iarg = iarg + 1
99 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
100 SELECT CASE (arg)
101 CASE ("-abc")
102 DO i = 1, 3
103 iarg = iarg + 1
104 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
105 READ (unit=arg, fmt=*, err=100) h(i, i)
106 END DO
107 have_cell_info = .true.
108 cycle arg_loop
109100 CALL abort_program(routinen, "Reading -abc arguments (3 reals are expected)")
110 CASE ("-cell")
111 DO i = 1, 3
112 DO j = 1, 3
113 iarg = iarg + 1
114 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
115 READ (unit=arg, fmt=*, err=101) h(j, i)
116 END DO
117 END DO
118 have_cell_info = .true.
119 cycle arg_loop
120101 CALL abort_program(routinen, "Reading -cell arguments (9 reals are expected)")
121 CASE ("-cell_file", "-cf")
122 iarg = iarg + 1
123 CALL get_command_argument(number=iarg, VALUE=cell_file_name, status=istat)
124 have_cell_file = .true.
125 have_cell_info = .true.
126 cycle arg_loop
127 CASE ("-df", "-dcd_file")
128 iarg = iarg + 1
129 CALL get_command_argument(number=iarg, VALUE=dcd_file_name, status=istat)
130 cycle arg_loop
131 CASE ("-debug", "-d")
132 debug = .true.
133 info = .true.
134 cycle arg_loop
135 CASE ("-displacements", "-disp")
136 print_atomic_displacements = .true.
137 cycle arg_loop
138 CASE ("-eo")
139 error_unit = output_unit
140 cycle arg_loop
141 CASE ("-first_frame", "-first", "-ff")
142 iarg = iarg + 1
143 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
144 READ (unit=arg, fmt=*, err=102) first_frame
145 IF (first_frame <= 0) THEN
146 CALL abort_program(routinen, "Invalid number for first frame specified: "// &
147 "first_frame must be greater than zero")
148 END IF
149 cycle arg_loop
150102 CALL abort_program(routinen, "Invalid number for first frame specified "// &
151 "(an integer number greater than zero is expected)")
152 CASE ("-help", "-h")
153 CALL print_help()
154 stop
155 CASE ("-info", "-i")
156 info = .true.
157 cycle arg_loop
158 CASE ("-last_frame", "-last", "-lf")
159 iarg = iarg + 1
160 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
161 READ (unit=arg, fmt=*, err=103) last_frame
162 IF (last_frame <= 0) THEN
163 CALL abort_program(routinen, "Invalid number for last frame specified: "// &
164 "last_frame must be greater than zero")
165 END IF
166 cycle arg_loop
167103 CALL abort_program(routinen, "Invalid number for last frame specified "// &
168 "(an integer number greater than zero is expected)")
169 CASE ("-pbc")
170 apply_pbc = .true.
171 pbc0 = .false.
172 cycle arg_loop
173 CASE ("-pbc0")
174 apply_pbc = .true.
175 pbc0 = .true.
176 cycle arg_loop
177 CASE ("-scaled_coordinates", "-sc")
178 print_scaled_coordinates = .true.
179 cycle arg_loop
180 CASE ("-scaled_pbc_coordinates", "-spc")
181 print_scaled_pbc_coordinates = .true.
182 cycle arg_loop
183 CASE ("-stride")
184 iarg = iarg + 1
185 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
186 READ (unit=arg, fmt=*, err=104) stride
187 IF (stride < 1) THEN
188 CALL abort_program(routinen, "Invalid stride for frame dump specified: stride must be greater than zero")
189 END IF
190 cycle arg_loop
191104 CALL abort_program(routinen, "Invalid stride for frame dump specified "// &
192 "(an integer number greater than 0 is expected)")
193 CASE ("-trace_atoms")
194 iarg = iarg + 1
195 CALL get_command_argument(number=iarg, VALUE=arg, status=istat)
196 READ (unit=arg, fmt=*, err=105) eps_out_of_box
197 IF (eps_out_of_box <= 0.0_dp) THEN
198 CALL abort_program(routinen, "Invalid threshold value for -trace_atoms flag specified")
199 END IF
200 trace_atoms = .true.
201 cycle arg_loop
202105 CALL abort_program(routinen, "Invalid threshold value for -trace_atoms flag specified")
203 CASE DEFAULT
204 IF (arg(1:1) == "-") THEN
205 CALL print_help()
206 CALL abort_program(routinen, "Unknown command line flag """//trim(arg)//""" found")
207 END IF
208 xyz_file_name = arg
209 EXIT arg_loop
210 END SELECT
211 END DO arg_loop
212
213 ! Check flag compatibility
214 IF (.NOT. have_cell_info) THEN
215 CALL abort_program(routinen, "No cell information available. Neither -abc, -cell, nor -cell_file flag found")
216 END IF
217 IF (first_frame > last_frame) THEN
218 CALL abort_program(routinen, "Number of first frame greater than number of last frame")
219 END IF
220 IF (.NOT. apply_pbc .AND. trace_atoms) THEN
221 CALL abort_program(routinen, "The -trace_atoms flag requires the specification of a -pbc flag")
222 END IF
223 IF (print_scaled_coordinates .AND. print_scaled_pbc_coordinates) THEN
224 CALL abort_program(routinen, "The -sc flag and the -spc flag are incompatible")
225 END IF
226 IF (.NOT. apply_pbc .AND. print_scaled_coordinates) THEN
227 CALL abort_program(routinen, "The -sc flag requires the specification of a -pbc flag")
228 END IF
229 IF (.NOT. apply_pbc .AND. print_scaled_pbc_coordinates) THEN
230 CALL abort_program(routinen, "The -spc flag requires the specification of a -pbc flag")
231 END IF
232
233 ! Open cell input file (if specified)
234 IF (have_cell_file) THEN
235 INQUIRE (file=cell_file_name, exist=exists)
236 IF (.NOT. exists) CALL abort_program(routinen, "The specified cell file <"// &
237 trim(cell_file_name)//"> does not exist")
238 OPEN (unit=cell_file_unit, &
239 file=cell_file_name, &
240 status="OLD", &
241 access="SEQUENTIAL", &
242 form="FORMATTED", &
243 position="REWIND", &
244 action="READ", &
245 iostat=istat)
246 IF (istat /= 0) CALL abort_program(routinen, "The cell file <"// &
247 trim(cell_file_name)//"> could not be opened")
248 IF (info) WRITE (unit=error_unit, fmt="(A)") "# Reading cell file : "//trim(cell_file_name)
249 END IF
250
251 ! Open XYZ input file
252 INQUIRE (file=xyz_file_name, exist=exists)
253 IF (.NOT. exists) CALL abort_program(routinen, "The specified XYZ file <"// &
254 trim(xyz_file_name)//"> does not exist")
255 OPEN (unit=xyz_file_unit, &
256 file=xyz_file_name, &
257 status="OLD", &
258 access="SEQUENTIAL", &
259 form="FORMATTED", &
260 position="REWIND", &
261 action="READ", &
262 iostat=istat)
263 IF (istat /= 0) CALL abort_program(routinen, "The XYZ file <"// &
264 trim(xyz_file_name)//"> could not be opened")
265 IF (info) WRITE (unit=error_unit, fmt="(A)") "# Reading XYZ file : "//trim(xyz_file_name)
266
267 ! Read the first two lines of the XYZ file
268 READ (unit=xyz_file_unit, fmt="(A)", iostat=istat) arg
269 IF (istat /= 0) THEN
270 CALL abort_program(routinen, "Reading the first line of the current frame from the XYZ file <"// &
271 trim(xyz_file_name)//"> failed")
272 END IF
273 IF (arg(1:1) == "#") THEN
274 READ (unit=arg, fmt=*) string, natom
275 ELSE
276 READ (unit=arg, fmt=*) natom
277 END IF
278 IF (istat /= 0) THEN
279 CALL abort_program(routinen, "Reading the number of atoms from the XYZ file failed")
280 END IF
281 READ (unit=xyz_file_unit, fmt="(A)", iostat=istat) remark_xyz
282 IF (istat /= 0) CALL abort_program(routinen, "Reading the second line from the XYZ file <"// &
283 trim(xyz_file_name)//"> failed")
284 rewind(unit=xyz_file_unit)
285
286 ! Open DCD output file
287 IF (len_trim(dcd_file_name) == 0) THEN
288 i = len_trim(xyz_file_name)
289 IF (xyz_file_name(i - 2:i) == "xyz") THEN
290 dcd_file_name = xyz_file_name(1:i - 3)//"dcd"
291 ELSE
292 dcd_file_name = xyz_file_name(1:i)//".dcd"
293 END IF
294 END IF
295 INQUIRE (file=dcd_file_name, exist=exists)
296 IF (exists) CALL abort_program(routinen, "The DCD file "// &
297 trim(dcd_file_name)//" exists already")
298 OPEN (unit=dcd_file_unit, &
299 file=dcd_file_name, &
300 status="UNKNOWN", &
301 access="SEQUENTIAL", &
302 form="UNFORMATTED", &
303 action="WRITE", &
304 iostat=istat)
305 IF (istat /= 0) CALL abort_program(routinen, "The unformatted DCD output file "// &
306 trim(dcd_file_name)//" could not be opened")
307 IF (info) WRITE (unit=error_unit, fmt="(A)") "# Writing DCD file : "//trim(dcd_file_name)
308
309 istep = 0
310 iskip = 0
311 dt = 0.0_dp
312
313 ! Write DCD file header
314 WRITE (unit=dcd_file_unit) "CORD", 0, istep, iskip, 0, 0, 0, 0, 0, 0, real(dt, kind=sp), &
315 1, 0, 0, 0, 0, 0, 0, 0, 0, 24
316 remark1 = "REMARK CORD"//" DCD file created by "//trim(version_info)
317 remark2 = "REMARK "//trim(adjustl(remark_xyz))
318 WRITE (unit=dcd_file_unit) 2, remark1(1:80), remark2(1:80)
319 WRITE (unit=dcd_file_unit) natom
320
321 ! Allocate work arrays
322 ALLOCATE (r(natom, 3), stat=istat)
323 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array r failed")
324 r(:, :) = 0.0_dp
325
326 ALLOCATE (atomic_label(natom), stat=istat)
327 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the vector atomic_label failed")
328 atomic_label(:) = ""
329
330 ! Loop over all frames in the XYZ file
331 frame_loop: DO
332
333 nframe = nframe + 1
334
335 IF (nframe < first_frame) THEN
336 dump_frame = .false.
337 ELSE
338 IF (modulo(nframe - first_frame, stride) == 0) THEN
339 dump_frame = .true.
340 ELSE
341 dump_frame = .false.
342 END IF
343 END IF
344
345 ! Read unit cell information, if available
346 IF (have_cell_file) THEN
347 DO
348 READ (unit=cell_file_unit, fmt=*, iostat=istat) arg
349 IF (istat < 0) EXIT frame_loop
350 IF (istat /= 0) THEN
351 CALL abort_program(routinen, "Reading line from cell file")
352 END IF
353 IF (arg(1:1) == "#") THEN
354 cycle
355 ELSE
356 backspace(unit=cell_file_unit)
357 EXIT
358 END IF
359 END DO
360 READ (unit=cell_file_unit, fmt=*, iostat=istat) istep, tstep, ((h(j, i), j=1, 3), i=1, 3)
361 IF (istat /= 0) THEN
362 CALL abort_program(routinen, "Reading information from cell file")
363 END IF
364 END IF
365
366 ! Initialise or update cell info
367 IF (have_cell_file .OR. (nframe == 1)) THEN
368 a(1:3) = h(1:3, 1)
369 b(1:3) = h(1:3, 2)
370 c(1:3) = h(1:3, 3)
371 abc(1) = norm2(a(1:3))
372 abc(2) = norm2(b(1:3))
373 abc(3) = norm2(c(1:3))
374 alpha = angle(b(1:3), c(1:3))*degree
375 beta = angle(a(1:3), c(1:3))*degree
376 gamma = angle(a(1:3), b(1:3))*degree
377 END IF
378
379 ! Read first line of the current frame in the XYZ file
380 READ (unit=xyz_file_unit, fmt="(A)", iostat=istat) arg
381 IF (istat < 0) EXIT frame_loop
382 IF (istat /= 0) THEN
383 CALL abort_program(routinen, "Reading the first line of the current frame from the XYZ file failed")
384 END IF
385 IF (arg(1:1) == "#") THEN
386 READ (unit=arg, fmt=*) string, natom
387 ELSE
388 READ (unit=arg, fmt=*) natom
389 END IF
390 IF (istat /= 0) THEN
391 CALL abort_program(routinen, "Reading the number of atoms from the XYZ file failed")
392 END IF
393 IF (natom /= SIZE(r, 1)) THEN
394 CALL abort_program(routinen, "Number of atoms changed for the current frame")
395 END IF
396
397 ! Read second line of the current frame in the XYZ file
398 READ (unit=xyz_file_unit, fmt="(A)", iostat=istat) remark_xyz
399 IF (istat /= 0) CALL abort_program(routinen, "Reading the second line from the XYZ file failed")
400
401 ! Optionally print some info
402 IF (info .AND. dump_frame) THEN
403 WRITE (unit=error_unit, fmt="(A,/,A,I0)") &
404 "#", "# Frame number : ", nframe
405 WRITE (unit=error_unit, fmt="(A,/,(A,F12.6))") &
406 "#", &
407 "# a [Angstrom] : ", abc(1), &
408 "# b [Angstrom] : ", abc(2), &
409 "# c [Angstrom] : ", abc(3), &
410 "# alpha [degree] : ", alpha, &
411 "# beta [degree] : ", beta, &
412 "# gamma [degree] : ", gamma
413 END IF
414
415 IF (info) THEN
416 WRITE (unit=output_unit, fmt="(T2,I0)") natom
417 WRITE (unit=output_unit, fmt="(A)") trim(adjustl(remark_xyz))
418 END IF
419
420 ! Read in the atomic positions of the current frame from the XYZ file
421 iatom = 0
422 DO
423 READ (unit=xyz_file_unit, fmt=*) arg
424 IF (arg(1:1) == "#") THEN
425 cycle
426 ELSE
427 backspace(unit=xyz_file_unit)
428 END IF
429 iatom = iatom + 1
430 READ (unit=xyz_file_unit, fmt=*, iostat=istat) atomic_label(iatom), r(iatom, 1:3)
431 IF (istat /= 0) THEN
432 message = ""
433 WRITE (unit=message, fmt="(A,I0,A,I0,A)") &
434 "Reading line ", iatom + 2, " of the current frame from XYZ file (atom ", iatom, ") failed"
435 CALL abort_program(routinen, trim(message))
436 END IF
437 CALL uppercase(atomic_label(iatom) (1:1))
438 IF (len_trim(atomic_label(iatom)) > 1) CALL lowercase(atomic_label(iatom) (2:2))
439 atomic_label(iatom) = trim(adjustl(atomic_label(iatom)))
440 IF (iatom == natom) EXIT
441 END DO
442
443 ! Store the atomic positions of the first frame in the current XYZ input file for
444 ! the output of the atomic displacements
445 IF ((nframe == 1) .AND. print_atomic_displacements) THEN
446 IF (.NOT. ALLOCATED(r0)) THEN
447 ALLOCATE (r0(natom, 3), stat=istat)
448 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array r0 failed")
449 END IF
450 r0(:, :) = r(:, :)
451 IF (.NOT. ALLOCATED(atomic_displacement)) THEN
452 ALLOCATE (atomic_displacement(natom), stat=istat)
453 IF (istat /= 0) THEN
454 CALL abort_program(routinen, "Allocation of the vector atomic_displacement failed")
455 END IF
456 END IF
457 atomic_displacement(:) = 0.0_dp
458 END IF
459
460 IF (dump_frame) THEN
461 ! Apply periodic boundary conditions before dumping the coordinates if requested
462 IF (apply_pbc) THEN
463 IF (.NOT. ALLOCATED(r_pbc)) THEN
464 ALLOCATE (r_pbc(natom, 3), stat=istat)
465 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array r_pbc failed")
466 r_pbc(:, :) = 0.0_dp
467 END IF
468 IF (.NOT. ALLOCATED(s)) THEN
469 ALLOCATE (s(natom, 3), stat=istat)
470 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array s failed")
471 s(:, :) = 0.0_dp
472 END IF
473 IF (.NOT. ALLOCATED(s_pbc)) THEN
474 ALLOCATE (s_pbc(natom, 3), stat=istat)
475 IF (istat /= 0) CALL abort_program(routinen, "Allocation of the array s_pbc failed")
476 s_pbc(:, :) = 0.0_dp
477 END IF
478 CALL pbc(r, r_pbc, s, s_pbc, h, hinv, debug, info, pbc0)
479 CALL write_out_of_box_atoms(atomic_label, r, s, eps_out_of_box, h)
480 ! Overwrite input coordinate with the PBCed coordinates for printing
481 r(:, :) = r_pbc(:, :)
482 END IF ! apply_pbc
483 ! Calculate the atomic displacements with respect to the first frame,
484 ! i.e. set of atomic positions
485 IF (print_atomic_displacements) THEN
486 DO iatom = 1, natom
487 atomic_displacement(iatom) = sqrt((r(iatom, 1) - r0(iatom, 1))**2 + &
488 (r(iatom, 2) - r0(iatom, 2))**2 + &
489 (r(iatom, 3) - r0(iatom, 3))**2)
490 END DO
491 END IF
492 IF (info) THEN
493 IF (print_scaled_coordinates) THEN
494 DO iatom = 1, natom
495 WRITE (unit=output_unit, fmt="(A5,3(1X,F14.6))") adjustl(atomic_label(iatom)), s(iatom, 1:3)
496 END DO
497 ELSE IF (print_scaled_pbc_coordinates) THEN
498 DO iatom = 1, natom
499 WRITE (unit=output_unit, fmt="(A5,3(1X,F14.6))") adjustl(atomic_label(iatom)), s_pbc(iatom, 1:3)
500 END DO
501 ELSE
502 DO iatom = 1, natom
503 WRITE (unit=output_unit, fmt="(A5,3(1X,F14.6))") adjustl(atomic_label(iatom)), r(iatom, 1:3)
504 END DO
505 END IF
506 END IF
507 ! Dump the cell information in DCD format
508 WRITE (unit=dcd_file_unit) abc(1), gamma, abc(2), beta, alpha, abc(3)
509 ! Dump the atomic coordinates in DCD format
510 DO i = 1, 3
511 WRITE (unit=dcd_file_unit) real(r(1:natom, i), kind=sp)
512 END DO
513 nframe_read = nframe_read + 1
514 END IF ! dump_frame
515
516 ! Exit loop and stop processing, if the last (requested) frame was encountered
517 IF (nframe >= last_frame) EXIT frame_loop
518
519 END DO frame_loop
520
521 nframe = nframe - 1
522
523 ! Close files
524 IF (have_cell_file) CLOSE (unit=cell_file_unit)
525 CLOSE (unit=dcd_file_unit)
526 CLOSE (unit=xyz_file_unit)
527
528 IF (info) THEN
529 WRITE (unit=error_unit, fmt="(A,/,A,I0)") &
530 "#", &
531 "# Frames processed : ", nframe_read
532 WRITE (unit=error_unit, fmt="(A)") &
533 "#", &
534 "# Normal termination of "//trim(version_info)
535 END IF
536
537 ! Cleanup
538 IF (ALLOCATED(atomic_label)) DEALLOCATE (atomic_label)
539 IF (ALLOCATED(atomic_displacement)) DEALLOCATE (atomic_displacement)
540 IF (ALLOCATED(r)) DEALLOCATE (r)
541 IF (ALLOCATED(r0)) DEALLOCATE (r0)
542 IF (ALLOCATED(r_pbc)) DEALLOCATE (r_pbc)
543 IF (ALLOCATED(s)) DEALLOCATE (s)
544 IF (ALLOCATED(s_pbc)) DEALLOCATE (s_pbc)
545
546CONTAINS
547
548! **************************************************************************************************
549!> \brief ...
550!> \param routine ...
551!> \param message ...
552! **************************************************************************************************
553 SUBROUTINE abort_program(routine, message)
554 ! Abort the program after printing an error message to stderr
555
556 CHARACTER(LEN=*), INTENT(IN) :: routine, message
557
558 CHARACTER(LEN=2*default_string_length) :: error_message
559
560 error_message = "*** ERROR in "//trim(routine)//": "//trim(message)//" ***"
561 WRITE (unit=default_error_unit, fmt="(/,A,/)") trim(error_message)
562 stop "*** ABNORMAL PROGRAM TERMINATION of xyz2dcd v1.0 ***"
563
564 END SUBROUTINE abort_program
565
566! **************************************************************************************************
567!> \brief ...
568!> \param a ...
569!> \param b ...
570!> \return ...
571! **************************************************************************************************
572 PURE FUNCTION angle(a, b) RESULT(angle_ab)
573 ! Calculation of the angle between the vectors a and b. The angle is returned in radians.
574
575 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: a, b
576 REAL(kind=dp) :: angle_ab
577
578 REAL(kind=dp), PARAMETER :: eps_geo = 1.0e-6_dp
579
580 REAL(kind=dp) :: length_of_a, length_of_b
581 REAL(kind=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
582
583 length_of_a = norm2(a)
584 length_of_b = norm2(b)
585
586 IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
587 a_norm(:) = a(:)/length_of_a
588 b_norm(:) = b(:)/length_of_b
589 angle_ab = acos(min(max(dot_product(a_norm, b_norm), -1.0_dp), 1.0_dp))
590 ELSE
591 angle_ab = 0.0_dp
592 END IF
593
594 END FUNCTION angle
595
596! **************************************************************************************************
597!> \brief ...
598!> \param a ...
599!> \param b ...
600!> \param c ...
601!> \param alpha ...
602!> \param beta ...
603!> \param gamma ...
604!> \param h ...
605! **************************************************************************************************
606 SUBROUTINE build_h_matrix(a, b, c, alpha, beta, gamma, h)
607 ! Calculate the h matrix of a simulation cell given the cell edge lengths a, b,
608 ! and c in Angstrom and the angles alpha (<b,c)), beta (<(a,c)), and gamma (<(a,b))
609 ! in degree.
610
611 REAL(kind=dp), INTENT(IN) :: a, b, c, alpha, beta, gamma
612 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: h
613
614 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_h_matrix'
615
616 REAL(kind=dp) :: cosa, cosb, cosg, sing
617
618 cosa = cos(alpha/degree)
619 IF (abs(cosa) < epsilon(0.0_dp)) cosa = 0.0_dp
620
621 cosb = cos(beta/degree)
622 IF (abs(cosb) < epsilon(0.0_dp)) cosb = 0.0_dp
623
624 cosg = cos(gamma/degree)
625 IF (abs(cosg) < epsilon(0.0_dp)) cosg = 0.0_dp
626
627 sing = sin(gamma/degree)
628 IF (abs(sing) < epsilon(0.0_dp)) sing = 0.0_dp
629
630 h(1, 1) = 1.0_dp
631 h(2, 1) = 0.0_dp
632 h(3, 1) = 0.0_dp
633
634 h(1, 2) = cosg
635 h(2, 2) = sing
636 h(3, 2) = 0.0_dp
637
638 h(1, 3) = cosb
639 h(2, 3) = (cosa - cosg*cosb)/sing
640 IF ((1.0_dp - h(1, 3)**2 - h(2, 3)**2) < 0.0_dp) THEN
641 CALL abort_program(routinen, "Build of the h matrix failed, check cell information")
642 END IF
643 h(3, 3) = sqrt(1.0_dp - h(1, 3)**2 - h(2, 3)**2)
644
645 h(:, 1) = a*h(:, 1)
646 h(:, 2) = b*h(:, 2)
647 h(:, 3) = c*h(:, 3)
648
649 END SUBROUTINE build_h_matrix
650
651! **************************************************************************************************
652!> \brief ...
653!> \param a ...
654!> \return ...
655! **************************************************************************************************
656 FUNCTION det_3x3(a) RESULT(det_a)
657 ! Returns the determinante of the 3x3 matrix a.
658
659 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: a
660 REAL(kind=dp) :: det_a
661
662 det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
663 a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
664 a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
665
666 END FUNCTION det_3x3
667
668! **************************************************************************************************
669!> \brief ...
670!> \param h ...
671!> \param hinv ...
672!> \param deth ...
673! **************************************************************************************************
674 SUBROUTINE invert_matrix_3x3(h, hinv, deth)
675 ! Calculate the inverse hinv and the determinant deth of the 3x3 matrix h.
676
677 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: h
678 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: hinv
679 REAL(kind=dp), INTENT(OUT) :: deth
680
681 CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_matrix_3x3'
682
683 deth = det_3x3(h)
684
685 ! Numerics
686 deth = abs(deth)
687 IF (deth < 1.0e-10_dp) THEN
688 CALL abort_program(routinen, "Invalid h matrix for cell found; det(h) < 1.0E-10")
689 END IF
690
691 hinv(1, 1) = (h(2, 2)*h(3, 3) - h(3, 2)*h(2, 3))/deth
692 hinv(2, 1) = (h(2, 3)*h(3, 1) - h(3, 3)*h(2, 1))/deth
693 hinv(3, 1) = (h(2, 1)*h(3, 2) - h(3, 1)*h(2, 2))/deth
694
695 hinv(1, 2) = (h(1, 3)*h(3, 2) - h(3, 3)*h(1, 2))/deth
696 hinv(2, 2) = (h(1, 1)*h(3, 3) - h(3, 1)*h(1, 3))/deth
697 hinv(3, 2) = (h(1, 2)*h(3, 1) - h(3, 2)*h(1, 1))/deth
698
699 hinv(1, 3) = (h(1, 2)*h(2, 3) - h(2, 2)*h(1, 3))/deth
700 hinv(2, 3) = (h(1, 3)*h(2, 1) - h(2, 3)*h(1, 1))/deth
701 hinv(3, 3) = (h(1, 1)*h(2, 2) - h(2, 1)*h(1, 2))/deth
702
703 END SUBROUTINE invert_matrix_3x3
704
705! **************************************************************************************************
706!> \brief ...
707!> \param string ...
708! **************************************************************************************************
709 SUBROUTINE lowercase(string)
710 ! Convert all letters in a string to lowercase
711 CHARACTER(LEN=*), INTENT(INOUT) :: string
712
713 INTEGER :: i, iascii
714
715 DO i = 1, len_trim(string)
716 iascii = ichar(string(i:i))
717 IF ((iascii >= 65) .AND. (iascii <= 90)) THEN
718 string(i:i) = char(iascii + 32)
719 END IF
720 END DO
721
722 END SUBROUTINE lowercase
723
724! **************************************************************************************************
725!> \brief ...
726!> \param r ...
727!> \param r_pbc ...
728!> \param s ...
729!> \param s_pbc ...
730!> \param h ...
731!> \param hinv ...
732!> \param debug ...
733!> \param info ...
734!> \param pbc0 ...
735! **************************************************************************************************
736 SUBROUTINE pbc(r, r_pbc, s, s_pbc, h, hinv, debug, info, pbc0)
737 ! Apply the periodic boundary conditions (PBC) to the Cartesian coordinate array
738 ! r given the cell edge lengths a, b, and c in Angstrom and the angles alpha (<b,c)),
739 ! beta (<(a,c)), and gamma (<(a,b)) in degree.
740 ! On output r_pbc is updated with the "PBCed" input coordinates, s with the scaled
741 ! input coordinates, and s_pbc with scaled "PBCed" coordinates.
742 ! If pbc0 is true then fold to the range [-l/2,+l/2[ (origin at box centre) else fold
743 ! to the range [0,l[ (origin at lower left box corner).
744
745 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: r
746 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: r_pbc, s, s_pbc
747 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: h
748 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: hinv
749 LOGICAL, INTENT(IN) :: debug, info, pbc0
750
751 CHARACTER(LEN=*), PARAMETER :: routinen = 'pbc'
752
753 INTEGER :: i, natom
754 LOGICAL :: orthorhombic
755 REAL(kind=dp) :: deth
756
757 natom = SIZE(r, 1)
758 IF (SIZE(r, 2) /= 3) CALL abort_program(routinen, "Array dimension for r must be 3")
759
760 orthorhombic = ((h(2, 1) == 0.0_dp) .AND. &
761 (h(3, 1) == 0.0_dp) .AND. &
762 (h(1, 2) == 0.0_dp) .AND. &
763 (h(3, 2) == 0.0_dp) .AND. &
764 (h(1, 3) == 0.0_dp) .AND. &
765 (h(2, 3) == 0.0_dp))
766
767 ! Build inverse of h
768 hinv(:, :) = 0.0_dp
769 CALL invert_matrix_3x3(h, hinv, deth)
770
771 IF (info) THEN
772 WRITE (unit=error_unit, fmt="(A)") "#"
773 IF (orthorhombic) THEN
774 WRITE (unit=error_unit, fmt="(A)") "# Cell symmetry : orthorhombic"
775 ELSE
776 WRITE (unit=error_unit, fmt="(A)") "# Cell symmetry : non-orthorhombic"
777 END IF
778 IF (debug) THEN
779 WRITE (unit=error_unit, fmt="(A)") "#"
780 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# / ", h(1, :), " \"
781 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# h = | ", h(2, :), " |"
782 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# \ ", h(3, :), " /"
783 WRITE (unit=error_unit, fmt="(A)") "#"
784 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# / ", hinv(1, :), " \"
785 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# Inv(h) = | ", hinv(2, :), " |"
786 WRITE (unit=error_unit, fmt="(A,3F12.6,A)") "# \ ", hinv(3, :), " /"
787 WRITE (unit=error_unit, fmt="(A)") "#"
788 WRITE (unit=error_unit, fmt="(A,F0.6)") "# det(h) = ", deth
789 END IF
790 END IF
791
792 ! Calculate scaled coordinates and wrap back all atoms, i.e. apply the PBC
793 IF (orthorhombic) THEN
794 ! Try to save some flops in the case of an orthorhombic box
795 DO i = 1, 3
796 s(:, i) = r(:, i)*hinv(i, i)
797 END DO
798 ELSE
799 s(:, :) = matmul(r(:, :), transpose(hinv(:, :)))
800 END IF
801
802 IF (pbc0) THEN
803 s_pbc(:, :) = s(:, :) - anint(s(:, :))
804 ELSE
805 s_pbc(:, :) = s(:, :) - floor(s(:, :))
806 END IF
807
808 IF (orthorhombic) THEN
809 DO i = 1, 3
810 r_pbc(:, i) = s_pbc(:, i)*h(i, i)
811 END DO
812 ELSE
813 r_pbc(:, :) = matmul(s_pbc(:, :), transpose(h(:, :)))
814 END IF
815
816 END SUBROUTINE pbc
817
818! **************************************************************************************************
819!> \brief ...
820! **************************************************************************************************
821 SUBROUTINE print_help()
822 ! Print the program flags for help
823
824 WRITE (unit=*, fmt="(T2,A)") &
825 "", &
826 "Program flags for "//trim(version_info)//":", &
827 "", &
828 " -abc <3 reals> : Cell vector lengths in Angstrom", &
829 " -cell <9 reals> : Cell vectors: a(x) a(y) a(z) b(x) b(y) b(z) c(x) c(y) c(z) in Angstrom", &
830 " -cell_file, -cf <file> : Name of the cell file in CP2K format", &
831 " -debug, -d : Print debug information", &
832 " -eo : Write standard output and standard error to the same logical unit", &
833 " -first_frame, -ff <int> : Number of the first frame which is dumped", &
834 " -help, -h : Print this information", &
835 " -info, -i : Print additional information for each frame (see also -debug flag)", &
836 " -last_frame, -lf <int> : Number of the last frame which is dumped", &
837 " -pbc : Apply the periodic boundary conditions (PBC) to each frame before it is dumped", &
838 " (origin at lower left)", &
839 " -pbc0 : Apply the periodic boundary conditions (PBC) to each frame before it is dumped", &
840 " (origin at box centre)", &
841 " -scaled_coordinates, -sc : Print the scaled coordinates", &
842 " -scaled_pbc_coordinates, -spc : Print the scaled coordinates after periodic boundary conditions (PBC) have been applied", &
843 " -stride <int> : Stride for frame dump (allows to skip frames, e.g. by dumping each 10th frame)", &
844 " -trace_atoms <real> : Print the atoms which left the simulation box given a threshold value in scaled units", &
845 ""
846
847 WRITE (unit=*, fmt="(T2,A)") &
848 "Usage examples:", &
849 "", &
850 " xyz2dcd <optional flags> <cell information: -abc <3 reals>, -cell <9 reals>, or -cell_file <file>> <XYZ file>", &
851 "", &
852 "Specific usage examples:", &
853 "", &
854 " xyz2dcd -abc 27.341 27.341 27.341 project-pos-1.xyz", &
855 " xyz2dcd -cell 27.341 0.0 0.0 0.0 27.341 0.0 0.0 0.0 27.341 project-pos-1.xyz", &
856 " xyz2dcd -cell_file project-1.cell project-pos-1.xyz", &
857 ""
858
859 WRITE (unit=*, fmt="(T2,A)") &
860 "Notes:", &
861 "", &
862 " - The -info and the -debug flags provide a more detailed output which is especially handy for tracing problems", &
863 " - The input coordinates and cell vectors should be in Angstrom", &
864 ""
865
866 END SUBROUTINE print_help
867
868! **************************************************************************************************
869!> \brief ...
870!> \param string ...
871! **************************************************************************************************
872 SUBROUTINE uppercase(string)
873 ! Convert all letters in a string to uppercase
874
875 CHARACTER(LEN=*), INTENT(INOUT) :: string
876
877 INTEGER :: i, iascii
878
879 DO i = 1, len_trim(string)
880 iascii = ichar(string(i:i))
881 IF ((iascii >= 97) .AND. (iascii <= 122)) THEN
882 string(i:i) = char(iascii - 32)
883 END IF
884 END DO
885
886 END SUBROUTINE uppercase
887
888! **************************************************************************************************
889!> \brief ...
890!> \param atomic_label ...
891!> \param r ...
892!> \param s ...
893!> \param eps_out_of_box ...
894!> \param h ...
895! **************************************************************************************************
896 SUBROUTINE write_out_of_box_atoms(atomic_label, r, s, eps_out_of_box, h)
897 ! Print a list of all atoms which have left the simulation box
898
899 CHARACTER(LEN=5), DIMENSION(:), INTENT(IN) :: atomic_label
900 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: r, s
901 REAL(kind=dp), INTENT(IN) :: eps_out_of_box
902 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: h
903
904 INTEGER :: i, iatom, natom, ncount
905 REAL(kind=dp) :: rl, s_max, s_min, sl
906 REAL(kind=dp), DIMENSION(3) :: dr, ds
907
908 ! Quick return, if no action is requested
909 IF (eps_out_of_box <= 0.0_dp) RETURN
910
911 s_max = 1.0_dp + eps_out_of_box
912 s_min = -eps_out_of_box
913 natom = SIZE(s, 1)
914 ncount = 0
915 DO iatom = 1, natom
916 IF (any(s(iatom, :) < s_min) .OR. &
917 any(s(iatom, :) > s_max)) THEN
918 ncount = ncount + 1
919 IF (ncount == 1) THEN
920 WRITE (unit=error_unit, fmt="(A)") &
921 "#", &
922 "# Atoms out of box:", &
923 "# Atom index label x y z |dr| |ds|"
924 END IF
925 ds(:) = s(iatom, :)
926 DO i = 1, 3
927 IF (s(iatom, i) < 0.0_dp) ds(i) = 0.0_dp
928 IF (s(iatom, i) >= 1.0_dp) ds(i) = 1.0_dp
929 END DO
930 ds(:) = s(iatom, :) - ds(:)
931 sl = sqrt(ds(1)**2 + ds(2)**2 + ds(3)**2)
932 dr(:) = matmul(h(:, :), ds(:))
933 rl = sqrt(dr(1)**2 + dr(2)**2 + dr(3)**2)
934 WRITE (unit=error_unit, fmt="(A,I10,1X,A5,5(1X,F14.6))") &
935 "# ", iatom, adjustr(atomic_label(iatom)), r(iatom, :), rl, sl
936 END IF
937 END DO
938 WRITE (unit=error_unit, fmt="(A,I0,A)") "# ", ncount, " atom(s) out of box"
939
940 END SUBROUTINE write_out_of_box_atoms
941
942END PROGRAM xyz2dcd
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
xyz2dcd
program xyz2dcd
Definition xyz2dcd.F:8