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