1PROGRAM dumpdcd 2 3! Copyright (C) 2012 - 2016 Matthias Krack (MK) 4! 5! Version: 3.0 6! Author: Matthias Krack (MK) 7! History: - Creation (13.02.2012,MK) 8! - XYZ file option added (13.02.2012,MK) 9! - Flags added for first and last frame and -o for output (24.05.2012,MK) 10! - XMOL flag added (25.05.2012,MK) 11! - PBC flag added (04.06.2012,MK) 12! - Stride flag added (05.06.2012,MK) 13! - Tracing of atoms added to detect the atoms which left the box (06.06.2012,MK) 14! - Keep input coordinates for further processing steps (15.06.2012,MK) 15! - VEL to CORD (-vel2cord flag) hack added (25.06.2012,MK) 16! - Added -displacement (-disp) flag (26.06.2012,MK) 17! - Dump the atomic displacements (CORD file) or temperatures (VEL file as x-coordinates of a DCD file (28.06.2012,MK) 18! - FRC to CORD (-frc2cord flag) hack added (28.01.2016,MK) 19! - Option -eformat added (15.02.2016,MK) 20! 21! Note: For -ekin a XYZ file is required to obtain the atomic labels. 22! The -info and the -debug flags provide a more detailed output which is especially handy for tracing problems. 23! The output in DCD format is in binary format. 24! The input coordinates should be in Angstrom. Velocities and forces are expected to be in atomic units. 25 26! Uncomment the following line if this module is available (e.g. with gfortran) and comment the corresponding variable declarations below 27! USE ISO_FORTRAN_ENV, ONLY: error_unit,input_unit,output_unit 28 29 IMPLICIT NONE 30 31 ! Comment the following lines if the ISO_FORTRAN_ENV is used (see above) 32 INTEGER, PARAMETER :: default_error_unit = 0,& 33 default_input_unit = 5,& 34 default_output_unit = 6 35 INTEGER :: error_unit = default_error_unit,& 36 input_unit = default_input_unit,& 37 output_unit = default_output_unit 38 ! End Comment 39 40 ! Parameters 41 CHARACTER(LEN=*), PARAMETER :: routine_name = "dumpdcd",& 42 version_info = routine_name//" v3.0 (31.01.2014, Matthias Krack)" 43 44 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14,200),& 45 sp = SELECTED_REAL_KIND(6,30) 46 INTEGER, PARAMETER :: default_string_length = 240,& 47 xyz_input_unit = 10 48 49 REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp 50 REAL(KIND=dp), PARAMETER :: angstrom = 0.52917720859_dp,& ! [a.u.] -> [Angstrom] 51 degree = 180.0_dp/pi,& ! [rad] -> [degree] 52 kelvin = 315774.647902944_dp,& ! [a.u.] -> [K] 53 massunit = 1822.88484264550_dp ! [u] -> [a.u.] 54 55 ! Variables 56 CHARACTER(LEN=4) :: id_dcd 57 CHARACTER(LEN=10) :: unit_string 58 CHARACTER(LEN=17) :: fmt_string 59 CHARACTER(LEN=default_string_length) :: arg,dcd_file_name,message,out_file_name,output_format,& 60 remark_xyz,string,xyz_file_name 61 CHARACTER(LEN=5), DIMENSION(:), ALLOCATABLE :: atomic_label 62 CHARACTER(LEN=80), DIMENSION(2) :: remark_dcd 63 INTEGER :: first_frame,have_unit_cell,i,iarg,& 64 iatom,iskip,istat,istep,last_frame,narg,& 65 natom_dcd,natom_xyz,ndcd_file,nframe,& 66 nframe_read,nremark,stride 67 LOGICAL :: apply_pbc,debug,dump_frame,eformat,ekin,eo,have_atomic_labels,& 68 info,opened,output_format_dcd,output_format_xmol,& 69 pbc0,print_atomic_displacements,print_scaled_coordinates,& 70 print_scaled_pbc_coordinates,trace_atoms,vel2cord 71 REAL(KIND=sp) :: dt 72 REAL(KIND=dp) :: a,alpha,b,beta,c,eps_out_of_box,gamma,tavg,tavg_frame 73 INTEGER, DIMENSION(16) :: idum 74 REAL(KIND=dp), DIMENSION(3) :: rdum 75 REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: atomic_displacement,atomic_mass,atomic_temperature 76 REAL(KIND=dp), DIMENSION(3,3) :: h,hinv 77 REAL(KIND=sp), DIMENSION(:,:), ALLOCATABLE :: r 78 REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: r_pbc,r0,s,s_pbc 79 80 apply_pbc = .FALSE. 81 pbc0 = .FALSE. 82 debug = .FALSE. 83 dump_frame = .TRUE. 84 eformat = .FALSE. 85 ekin = .FALSE. 86 eo = .FALSE. 87 info = .FALSE. 88 trace_atoms = .FALSE. 89 vel2cord = .FALSE. 90 print_atomic_displacements = .FALSE. 91 print_scaled_coordinates = .FALSE. 92 print_scaled_pbc_coordinates = .FALSE. 93 first_frame = 1 94 last_frame = 1000000 ! Hard limit of 1 Mio frames in total 95 stride = 1 96 ndcd_file = 0 97 nframe = 0 98 nframe_read = 0 99 nremark = 0 100 have_unit_cell = 0 101 have_atomic_labels = .FALSE. 102 idum(:) = 0 103 rdum(:) = 0.0_dp 104 id_dcd = "" 105 dcd_file_name = "" 106 out_file_name = "" 107 xyz_file_name = "" 108 output_format = "default" 109 output_format_dcd = .FALSE. 110 output_format_xmol = .FALSE. 111 remark_dcd(:) = "" 112 remark_xyz = "" 113 fmt_string = "" 114 dt = 0.0_sp 115 a = 0.0_dp 116 b = 0.0_dp 117 c = 0.0_dp 118 alpha = 0.0_dp 119 beta = 0.0_dp 120 gamma = 0.0_dp 121 eps_out_of_box = -HUGE(0.0_dp) 122 tavg = 0.0_dp 123 tavg_frame = 0.0_dp 124 125 narg = command_argument_count() 126 127 IF (narg == 0) THEN 128 CALL print_help() 129 CALL abort_program(routine_name,"No input file(s) specified") 130 END IF 131 132 iarg = 0 133 134 dcd_file_loop: DO 135 136 iarg = iarg + 1 137 138 CALL get_command_argument(NUMBER=iarg,VALUE=arg,STATUS=istat) 139 140 SELECT CASE (arg) 141 CASE ("-debug","-d") 142 debug = .TRUE. 143 info = .TRUE. 144 CYCLE dcd_file_loop 145 CASE ("-displacements","-disp") 146 print_atomic_displacements = .TRUE. 147 CYCLE dcd_file_loop 148 CASE ("-eformat") 149 eformat = .TRUE. 150 CYCLE dcd_file_loop 151 CASE ("-ekin") 152 ekin = .TRUE. 153 CYCLE dcd_file_loop 154 CASE ("-eo") 155 eo = .TRUE. 156 CYCLE dcd_file_loop 157 CASE ("-first_frame","-first","-ff") 158 iarg = iarg + 1 159 CALL get_command_argument(NUMBER=iarg,VALUE=arg,STATUS=istat) 160 READ (UNIT=arg,FMT=*,ERR=100) first_frame 161 IF (first_frame <= 0) THEN 162 CALL abort_program(routine_name,"Invalid number for first frame specified: "//& 163 "first_frame must be greater than zero") 164 END IF 165 CYCLE dcd_file_loop 166100 CALL abort_program(routine_name,"Invalid number for first frame specified "//& 167 "(an integer number greater than zero is expected)") 168 CASE ("-help","-h") 169 CALL print_help() 170 STOP 171 CASE ("-info","-i") 172 info = .TRUE. 173 CYCLE dcd_file_loop 174 CASE ("-last_frame","-last","-lf") 175 iarg = iarg + 1 176 CALL get_command_argument(NUMBER=iarg,VALUE=arg,STATUS=istat) 177 READ (UNIT=arg,FMT=*,ERR=101) last_frame 178 IF (last_frame <= 0) THEN 179 CALL abort_program(routine_name,"Invalid number for last frame specified: "//& 180 "last_frame must be greater than zero") 181 END IF 182 CYCLE dcd_file_loop 183101 CALL abort_program(routine_name,"Invalid number for last frame specified "//& 184 "(an integer number greater than zero is expected)") 185 CASE ("-o","-output") 186 iarg = iarg + 1 187 CALL get_command_argument(NUMBER=iarg,VALUE=out_file_name,STATUS=istat) 188 CYCLE dcd_file_loop 189 CASE ("-output_format","-of") 190 iarg = iarg + 1 191 CALL get_command_argument(NUMBER=iarg,VALUE=output_format,STATUS=istat) 192 CALL uppercase(output_format) 193 SELECT CASE (output_format) 194 CASE ("DCD") 195 output_format_dcd = .TRUE. 196 output_format_xmol = .FALSE. 197 CASE ("XMOL") 198 output_format_dcd = .FALSE. 199 output_format_xmol = .TRUE. 200 CASE DEFAULT 201 CALL abort_program(routine_name,"Invalid output format type specified") 202 END SELECT 203 CYCLE dcd_file_loop 204 CASE ("-pbc") 205 apply_pbc = .TRUE. 206 pbc0 = .FALSE. 207 CYCLE dcd_file_loop 208 CASE ("-pbc0") 209 apply_pbc = .TRUE. 210 pbc0 = .TRUE. 211 CYCLE dcd_file_loop 212 CASE ("-scaled_coordinates","-sc") 213 print_scaled_coordinates = .TRUE. 214 CYCLE dcd_file_loop 215 CASE ("-scaled_pbc_coordinates","-spc") 216 print_scaled_pbc_coordinates = .TRUE. 217 CYCLE dcd_file_loop 218 CASE ("-stride") 219 iarg = iarg + 1 220 CALL get_command_argument(NUMBER=iarg,VALUE=arg,STATUS=istat) 221 READ (UNIT=arg,FMT=*,ERR=104) stride 222 IF (stride < 1) THEN 223 CALL abort_program(routine_name,"Invalid stride for frame dump specified: stride must be greater than zero") 224 END IF 225 CYCLE dcd_file_loop 226104 CALL abort_program(routine_name,"Invalid stride for frame dump specified "//& 227 "(an integer number greater than 0 is expected)") 228 CASE ("-trace_atoms") 229 iarg = iarg + 1 230 CALL get_command_argument(NUMBER=iarg,VALUE=arg,STATUS=istat) 231 READ (UNIT=arg,FMT=*,ERR=108) eps_out_of_box 232 IF (eps_out_of_box <= 0.0_dp) THEN 233 CALL abort_program(routine_name,"Invalid threshold value for -trace_atoms flag specified") 234 END IF 235 trace_atoms = .TRUE. 236 CYCLE dcd_file_loop 237108 CALL abort_program(routine_name,"Invalid threshold value for -trace_atoms flag specified") 238 CASE ("-vel2cord","-v2c","-frc2cord","-f2c") 239 vel2cord = .TRUE. 240 CYCLE dcd_file_loop 241 CASE ("-xyz","xyz_file") 242 iarg = iarg + 1 243 CALL get_command_argument(NUMBER=iarg,VALUE=xyz_file_name,STATUS=istat) 244 have_atomic_labels = .TRUE. 245 CYCLE dcd_file_loop 246 CASE DEFAULT 247 IF (arg(1:1) == "-") THEN 248 CALL print_help() 249 CALL abort_program(routine_name,"Unknown command line flag """//TRIM(arg)//""" found") 250 END IF 251 dcd_file_name = arg 252 END SELECT 253 254 ! Check flag compatibility 255 IF (first_frame > last_frame) THEN 256 CALL abort_program(routine_name,"Number of first frame greater than number of last frame") 257 END IF 258 IF ((.NOT.have_atomic_labels).AND.output_format_xmol) THEN 259 CALL abort_program(routine_name,"The output format XMOL requires a valid xyz file (-xyz flag)") 260 END IF 261 IF (output_format_xmol.AND.ekin) THEN 262 CALL abort_program(routine_name,"Output format XMOL and the -ekin flag are incompatible") 263 END IF 264 IF (output_format_xmol.AND.print_atomic_displacements) THEN 265 CALL abort_program(routine_name,"Output format XMOL and the -displacements flag are incompatible") 266 END IF 267 IF (output_format_xmol.AND.eo) THEN 268 CALL abort_program(routine_name,"The -eo flag is incompatible with the output format XMOL") 269 END IF 270 IF (.NOT.have_atomic_labels.AND.ekin) THEN 271 CALL abort_program(routine_name,"ekin flag requires also the specification of a valid xyz file (-xyz flag)") 272 END IF 273 IF (.NOT.apply_pbc.AND.trace_atoms) THEN 274 CALL abort_program(routine_name,"The -trace_atoms flag requires the specification of a -pbc flag") 275 END IF 276 IF (ekin.AND.print_atomic_displacements) THEN 277 CALL abort_program(routine_name,"The -ekin flag and the -displacements flag are incompatible") 278 END IF 279 IF (print_scaled_coordinates.AND.print_scaled_pbc_coordinates) THEN 280 CALL abort_program(routine_name,"The -sc flag and the -spc flag are incompatible") 281 END IF 282 IF (.NOT.apply_pbc.AND.print_scaled_coordinates) THEN 283 CALL abort_program(routine_name,"The -sc flag requires the specification of a -pbc flag") 284 END IF 285 IF (.NOT.apply_pbc.AND.print_scaled_pbc_coordinates) THEN 286 CALL abort_program(routine_name,"The -spc flag requires the specification of a -pbc flag") 287 END IF 288 289 ! Use optionally an ES format 290 IF (eformat) THEN 291 fmt_string = "(A5,3(1X,ES14.6))" 292 ELSE 293 fmt_string = "(A5,3(1X, F14.6))" 294 END IF 295 296 ! Open output units as requested 297 IF (output_format_dcd) THEN 298 ! Set default output file name if no file name was specified 299 IF (LEN_TRIM(out_file_name) == 0) out_file_name = "output.dcd" 300 ! Check if a new output file name was specified, if yes then close the old unit 301 INQUIRE (UNIT=output_unit,NAME=string) 302 IF (TRIM(string) /= TRIM(out_file_name)) CLOSE (UNIT=output_unit) 303 INQUIRE (UNIT=output_unit,OPENED=opened) 304 IF (.NOT.opened) THEN 305 OPEN (UNIT=output_unit,& 306 FILE=out_file_name,& 307 STATUS="UNKNOWN",& 308 ACCESS="SEQUENTIAL",& 309 FORM="UNFORMATTED",& 310 ACTION="WRITE",& 311 IOSTAT=istat) 312 IF (istat /= 0) CALL abort_program(routine_name,"The unformatted output file could not be opened") 313 END IF 314 ELSE 315 IF (eo) error_unit = output_unit 316 IF (LEN_TRIM(out_file_name) > 0) THEN 317 ! Check if a new output file name was specified, if yes then close the old unit 318 INQUIRE (UNIT=output_unit,NAME=string) 319 IF (TRIM(string) /= TRIM(out_file_name)) CLOSE (UNIT=output_unit) 320 INQUIRE (UNIT=output_unit,OPENED=opened) 321 IF (.NOT.opened) THEN 322 OPEN (UNIT=output_unit,& 323 FILE=out_file_name,& 324 STATUS="UNKNOWN",& 325 ACCESS="SEQUENTIAL",& 326 FORM="FORMATTED",& 327 POSITION="REWIND",& 328 ACTION="WRITE",& 329 IOSTAT=istat) 330 IF (istat /= 0) CALL abort_program(routine_name,"The formatted output file could not be opened") 331 END IF 332 END IF 333 END IF 334 335 ! Avoid reading from and writing to the same file 336 IF (TRIM(dcd_file_name) == TRIM(out_file_name)) THEN 337 CALL abort_program(routine_name,"Input and output file name cannot be the same") 338 END IF 339 340 ! Read reference xyz input file if requested 341 IF (have_atomic_labels) THEN 342 string = "" 343 ! Check if a new XYZ input file name was specified, if yes then close the old unit 344 INQUIRE (UNIT=xyz_input_unit,NAME=string) 345 IF (TRIM(string) /= TRIM(xyz_file_name)) CLOSE (UNIT=xyz_input_unit) 346 INQUIRE (UNIT=xyz_input_unit,OPENED=opened) 347 ! Read new XYZ input file if needed and update reference information 348 IF (.NOT.opened) THEN 349 OPEN (UNIT=xyz_input_unit,& 350 FILE=xyz_file_name,& 351 STATUS="OLD",& 352 ACCESS="SEQUENTIAL",& 353 FORM="FORMATTED",& 354 POSITION="REWIND",& 355 ACTION="READ",& 356 IOSTAT=istat) 357 IF (istat /= 0) CALL abort_program(routine_name,"The XYZ file could not be opened") 358 IF (info) WRITE (UNIT=error_unit,FMT="(A)") "#","# Reading XYZ file: "//TRIM(xyz_file_name) 359 READ (UNIT=xyz_input_unit,FMT="(A)",IOSTAT=istat) arg 360 IF (istat /= 0) THEN 361 CALL abort_program(routine_name,"Reading line 1 of the XYZ file (number of atoms) failed") 362 END IF 363 IF (arg(1:1) == "#") THEN 364 READ (UNIT=arg,FMT=*,IOSTAT=istat) string,natom_xyz 365 ELSE 366 READ (UNIT=arg,FMT=*,IOSTAT=istat) natom_xyz 367 END IF 368 IF (istat /= 0) THEN 369 CALL abort_program(routine_name,"Reading line 1 of the XYZ file (number of atoms) failed") 370 END IF 371 IF (ALLOCATED(atomic_label)) THEN 372 DEALLOCATE (atomic_label,STAT=istat) 373 IF (istat /= 0) CALL abort_program(routine_name,"Deallocation of the vector atomic_label failed") 374 END IF 375 ALLOCATE (atomic_label(natom_xyz),STAT=istat) 376 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the vector atomic_label failed") 377 atomic_label(:) = "" 378 IF (ekin) THEN 379 IF (ALLOCATED(atomic_mass)) THEN 380 DEALLOCATE (atomic_mass,STAT=istat) 381 IF (istat /= 0) CALL abort_program(routine_name,"Deallocation of the vector atomic_mass failed") 382 END IF 383 ALLOCATE (atomic_mass(natom_xyz),STAT=istat) 384 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the vector atomic_mass failed") 385 END IF 386 READ (UNIT=xyz_input_unit,FMT="(A)",IOSTAT=istat) remark_xyz 387 IF (istat /= 0) CALL abort_program(routine_name,"Reading line 2 of the XYZ file (remark line) failed") 388 iatom = 0 389 DO 390 READ (UNIT=xyz_input_unit,FMT=*,IOSTAT=istat) arg 391 IF (arg(1:1) == "#") THEN 392 CYCLE 393 ELSE 394 BACKSPACE (UNIT=xyz_input_unit) 395 END IF 396 iatom = iatom + 1 397 READ (UNIT=xyz_input_unit,FMT=*,IOSTAT=istat) atomic_label(iatom),rdum(1:3) 398 IF (istat /= 0) THEN 399 message = "" 400 WRITE (UNIT=message,FMT="(A,I0,A,I0,A)")& 401 "Reading line ",iatom+2," of the XYZ file (atom ",iatom,") failed" 402 CALL abort_program(routine_name,TRIM(message)) 403 END IF 404 CALL uppercase(atomic_label(iatom)) 405 IF (LEN_TRIM(atomic_label(iatom)) > 1) CALL lowercase(atomic_label(iatom)(2:2)) 406 atomic_label(iatom) = TRIM(ADJUSTL(atomic_label(iatom))) 407 IF (ekin) atomic_mass(iatom) = get_atomic_mass(atomic_label(iatom)) 408 IF (iatom == natom_xyz) EXIT 409 END DO 410 IF (info) THEN 411 WRITE (UNIT=error_unit,FMT="(A,I0)")& 412 "# Number of atoms : ",natom_xyz,& 413 "# Remark : "//TRIM(ADJUSTL(remark_xyz)) 414 END IF 415 END IF 416 END IF 417 418 ! Increment DCD file counter 419 ndcd_file = ndcd_file + 1 420 421 ! Open input file in DCD format 422 OPEN (UNIT=input_unit,& 423 FILE=dcd_file_name,& 424 STATUS="OLD",& 425 FORM="UNFORMATTED",& 426 POSITION="REWIND",& 427 ACTION="READ",& 428 IOSTAT=istat) 429 IF (istat /= 0) CALL abort_program(routine_name,"The DCD file could not be opened") 430 IF (info) THEN 431 WRITE (UNIT=error_unit,FMT="(A,/,(A,I0))")& 432 "#",& 433 "# DCD file number : ",ndcd_file,& 434 "# Reading DCD file: "//TRIM(dcd_file_name) 435 END IF 436 437 ! Read 1st record of DCD file 438 READ (UNIT=input_unit) id_dcd,idum(1),istep,iskip,idum(2:7),dt,have_unit_cell,idum(8:16) 439 IF (info) THEN 440 WRITE (UNIT=error_unit,FMT="(A,2(/,A,I0),/,A,F9.3,/,A,I0)")& 441 "# DCD id string : "//id_dcd,& 442 "# Step : ",istep,& 443 "# Print frequency : ",iskip,& 444 "# Time step [fs] : ",dt,& 445 "# Unit cell : ",have_unit_cell 446 END IF 447 448 IF (ekin.AND.TRIM(ADJUSTL(id_dcd)) /= "VEL") THEN 449 CALL abort_program(routine_name,"ekin flag requires a DCD file with VELocities") 450 END IF 451 IF (apply_pbc.AND.(have_unit_cell /= 1)) THEN 452 CALL abort_program(routine_name,"pbc flags require that unit cell information is available") 453 END IF 454 455 IF ((TRIM(ADJUSTL(id_dcd)) == "FRC").OR.(TRIM(ADJUSTL(id_dcd)) == "VEL")) THEN 456 unit_string = "[a.u.]" 457 IF (apply_pbc) CALL abort_program(routine_name,"pbc flags require a DCD file with COoRDinates") 458 ELSE IF (TRIM(ADJUSTL(id_dcd)) == "CORD") THEN 459 unit_string = "[Angstrom]" 460 ELSE 461 CALL abort_program(routine_name,"Unknown DCD id found (use -debug or -info flag for details)") 462 END IF 463 464 ! Read 2nd record of DCD file 465 READ (UNIT=input_unit) nremark,remark_dcd(1),remark_dcd(2) 466 IF (info) THEN 467 DO i=1,nremark 468 WRITE (UNIT=error_unit,FMT="(A,I1,A)")& 469 "# Remark ",i," : "//TRIM(remark_dcd(i)) 470 END DO 471 END IF 472 473 ! Read 3rd record of DCD file 474 READ (UNIT=input_unit) natom_dcd 475 IF (info) THEN 476 WRITE (UNIT=error_unit,FMT="(A,I0)")& 477 "# Number of atoms : ",natom_dcd 478 END IF 479 480 IF (have_atomic_labels) THEN 481 IF (natom_dcd /= natom_xyz) CALL abort_program(routine_name,"Number of atoms in XYZ and DCD file differ") 482 ELSE 483 IF (.NOT.ALLOCATED(atomic_label)) THEN 484 ALLOCATE (atomic_label(natom_dcd),STAT=istat) 485 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the vector atomic_label failed") 486 atomic_label(:) = "" 487 END IF 488 END IF 489 490 ! Dump output in DCD format => dcd2dcd functionality 491 IF (output_format_dcd) THEN 492 IF (vel2cord) id_dcd = "CORD" 493 ! Note, dt is REAL*4 and the rest is INTEGER*4 494 WRITE (UNIT=output_unit) id_dcd,idum(1),istep,iskip,idum(2:7),dt,have_unit_cell,idum(8:16) 495 WRITE (UNIT=output_unit) nremark,remark_dcd(1),remark_dcd(2) 496 WRITE (UNIT=output_unit) natom_dcd 497 END IF 498 499 frame_loop: DO 500 501 nframe = nframe + 1 502 503 IF (nframe < first_frame) THEN 504 dump_frame = .FALSE. 505 ELSE 506 IF (MODULO(nframe-first_frame,stride) == 0) THEN 507 dump_frame = .TRUE. 508 ELSE 509 dump_frame = .FALSE. 510 END IF 511 END IF 512 513 ! Read unit cell information, if available 514 IF (have_unit_cell == 1) THEN 515 READ (UNIT=input_unit,IOSTAT=istat) a,gamma,b,beta,alpha,c 516 IF (istat < 0) EXIT frame_loop 517 END IF 518 519 IF ((info.OR.trace_atoms).AND.dump_frame) THEN 520 WRITE (UNIT=error_unit,FMT="(A,/,A,I0)")& 521 "#","# Frame number : ",nframe 522 END IF 523 524 ! Print unit cell information, if available 525 IF (info.AND.dump_frame) THEN 526 IF (have_unit_cell == 1) THEN 527 WRITE (UNIT=error_unit,FMT="(A,/,(A,F12.6))")& 528 "#",& 529 "# a [Angstrom] : ",a,& 530 "# b [Angstrom] : ",b,& 531 "# c [Angstrom] : ",c,& 532 "# alpha [degree] : ",alpha,& 533 "# beta [degree] : ",beta,& 534 "# gamma [degree] : ",gamma 535 END IF 536 END IF 537 538 ! Allocate the array for the current atomic positions if needed 539 IF (.NOT.ALLOCATED(r)) THEN 540 ALLOCATE (r(natom_dcd,3),STAT=istat) 541 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the array r failed") 542 END IF 543 544 ! Read in the atomic positions of the current frame 545 READ (UNIT=input_unit,IOSTAT=istat) r(1:natom_dcd,1) 546 IF (istat < 0) EXIT frame_loop 547 READ (UNIT=input_unit) r(1:natom_dcd,2) 548 READ (UNIT=input_unit) r(1:natom_dcd,3) 549 550 ! Store the atomic positions of the first frame in the current DCD input file for 551 ! the output of the atomic displacements 552 IF ((nframe == 1).AND.print_atomic_displacements) THEN 553 IF (.NOT.ALLOCATED(r0)) THEN 554 ALLOCATE (r0(natom_dcd,3),STAT=istat) 555 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the array r0 failed") 556 END IF 557 r0(:,:) = r(:,:) 558 IF (.NOT.ALLOCATED(atomic_displacement)) THEN 559 ALLOCATE (atomic_displacement(natom_dcd),STAT=istat) 560 IF (istat /= 0) THEN 561 CALL abort_program(routine_name,"Allocation of the vector atomic_displacement failed") 562 END IF 563 END IF 564 atomic_displacement(:) = 0.0_dp 565 END IF 566 567 IF (ekin.AND.TRIM(ADJUSTL(id_dcd)) == "VEL") THEN 568 569 ! Dump the kinetic energy of each as an "atomic temperature" 570 IF (dump_frame) THEN 571 572 IF (.NOT.ALLOCATED(atomic_temperature)) THEN 573 ALLOCATE (atomic_temperature(natom_dcd),STAT=istat) 574 IF (istat /= 0) THEN 575 CALL abort_program(routine_name,"Allocation of the vector atomic_temperature failed") 576 END IF 577 atomic_temperature(:) = 0.0_dp 578 END IF 579 580 IF (info) THEN 581 WRITE (UNIT=error_unit,FMT="(A)")& 582 "#",& 583 "# Temperature [K]" 584 END IF 585 586 tavg_frame = 0.0_dp 587 DO iatom=1,natom_dcd 588 atomic_temperature(iatom) = atomic_mass(iatom)*(r(iatom,1)*r(iatom,1) +& 589 r(iatom,2)*r(iatom,2) +& 590 r(iatom,3)*r(iatom,3))*kelvin/3.0_dp 591 tavg_frame = tavg_frame + atomic_temperature(iatom) 592 END DO 593 tavg_frame = tavg_frame/REAL(natom_dcd,KIND=dp) 594 595 IF (output_format_dcd) THEN 596 IF (have_unit_cell == 1) THEN 597 WRITE (UNIT=output_unit) a,gamma,b,beta,alpha,c 598 END IF 599 ! This is a hack for VMD: dump the atomic temperatures as the x-coordinates 600 ! of a DCD VEL file 601 WRITE (UNIT=output_unit) REAL(atomic_temperature(:),KIND=sp) 602 ! The y and z coordinates are filled with zeros 603 atomic_temperature(:) = 0.0_dp 604 WRITE (UNIT=output_unit) REAL(atomic_temperature(:),KIND=sp) 605 WRITE (UNIT=output_unit) REAL(atomic_temperature(:),KIND=sp) 606 ELSE 607 DO iatom=1,natom_dcd 608 WRITE (UNIT=output_unit,FMT="(A5,5X,F25.3)") ADJUSTL(atomic_label(iatom)),atomic_temperature 609 END DO 610 END IF 611 612 IF (info) THEN 613 WRITE (UNIT=error_unit,FMT="(A,F12.3)")& 614 "# T [K] this frame: ",tavg_frame 615 END IF 616 617 tavg = tavg + tavg_frame 618 619 END IF ! dump_frame 620 621 ELSE 622 623 IF (dump_frame) THEN 624 625 ! Apply periodic boundary conditions before dumping the coordinates if requested 626 IF (apply_pbc) THEN 627 IF (.NOT.ALLOCATED(r_pbc)) THEN 628 ALLOCATE (r_pbc(natom_dcd,3),STAT=istat) 629 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the array r_pbc failed") 630 END IF 631 IF (.NOT.ALLOCATED(s)) THEN 632 ALLOCATE (s(natom_dcd,3),STAT=istat) 633 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the array s failed") 634 END IF 635 IF (.NOT.ALLOCATED(s_pbc)) THEN 636 ALLOCATE (s_pbc(natom_dcd,3),STAT=istat) 637 IF (istat /= 0) CALL abort_program(routine_name,"Allocation of the array s_pbc failed") 638 END IF 639 CALL pbc(r,r_pbc,s,s_pbc,a,b,c,alpha,beta,gamma,debug,info,pbc0,h,hinv) 640 CALL write_out_of_box_atoms(atomic_label,r,s,eps_out_of_box,h) 641 ! Overwrite input coordinate with the PBCed coordinates for printing 642 r(:,:) = REAL(r_pbc(:,:),KIND=sp) 643 END IF ! apply_pbc 644 645 ! Calculate the atomic displacements with respect to the first frame, 646 ! i.e. set of atomic positions 647 IF (print_atomic_displacements) THEN 648 DO iatom=1,natom_dcd 649 atomic_displacement(iatom) = SQRT((r(iatom,1) - r0(iatom,1))**2 +& 650 (r(iatom,2) - r0(iatom,2))**2 +& 651 (r(iatom,3) - r0(iatom,3))**2) 652 END DO 653 END IF 654 655 ! Dump the atomic coordinates in DCD or XMOL/MOLDEN format 656 IF (output_format_dcd) THEN 657 IF (have_unit_cell == 1) THEN 658 WRITE (UNIT=output_unit) a,gamma,b,beta,alpha,c 659 END IF 660 IF (print_atomic_displacements) THEN 661 ! This is a hack for VMD: dump the atomic displacements as the x-coordinates 662 ! of a DCD CORD file 663 WRITE (UNIT=output_unit) REAL(atomic_displacement(:),KIND=sp) 664 ! The y and z coordinates are filled with zeros 665 atomic_displacement(:) = 0.0_dp 666 WRITE (UNIT=output_unit) REAL(atomic_displacement(:),KIND=sp) 667 WRITE (UNIT=output_unit) REAL(atomic_displacement(:),KIND=sp) 668 ELSE 669 ! Dump DCD file information 670 DO i=1,3 671 WRITE (UNIT=output_unit) r(:,i) 672 END DO 673 END IF 674 ELSE 675 IF (print_atomic_displacements) THEN 676 IF (info) THEN 677 WRITE (UNIT=error_unit,FMT="(A,2(/,A),T15,A)")& 678 "#",& 679 "# Displacements :",& 680 "# "//TRIM(unit_string),"|r|" 681 END IF 682 IF (have_atomic_labels) THEN 683 DO iatom=1,natom_dcd 684 WRITE (UNIT=output_unit,FMT="(A5,5X,F25.6)") ADJUSTL(atomic_label(iatom)),atomic_displacement(iatom) 685 END DO 686 ELSE 687 DO iatom=1,natom_dcd 688 WRITE (UNIT=output_unit,FMT="(10X,F25.6)") atomic_displacement(iatom) 689 END DO 690 END IF 691 ELSE 692 IF (output_format_xmol) THEN 693 WRITE (UNIT=output_unit,FMT="(T2,I0,/,A,2(I0,A))")& 694 natom_dcd,"Frame: ",nframe_read + 1,", Step: ",nframe,", "//TRIM(remark_xyz) 695 DO iatom=1,natom_dcd 696 WRITE (UNIT=output_unit,FMT=fmt_string) ADJUSTL(atomic_label(iatom)),r(iatom,1:3) 697 END DO 698 ELSE 699 IF (info) THEN 700 WRITE (UNIT=error_unit,FMT="(A,T20,A)")& 701 "# "//TRIM(unit_string),"x y z" 702 END IF 703 IF (print_scaled_coordinates) THEN 704 DO iatom=1,natom_dcd 705 WRITE (UNIT=output_unit,FMT=fmt_string) ADJUSTL(atomic_label(iatom)),s(iatom,1:3) 706 END DO 707 ELSE IF (print_scaled_pbc_coordinates) THEN 708 DO iatom=1,natom_dcd 709 WRITE (UNIT=output_unit,FMT=fmt_string) ADJUSTL(atomic_label(iatom)),s_pbc(iatom,1:3) 710 END DO 711 ELSE 712 DO iatom=1,natom_dcd 713 WRITE (UNIT=output_unit,FMT=fmt_string) ADJUSTL(atomic_label(iatom)),r(iatom,1:3) 714 END DO 715 END IF 716 END IF 717 END IF 718 END IF ! output_format_dcd 719 720 END IF ! dump_frame 721 722 END IF 723 724 IF (dump_frame) nframe_read = nframe_read + 1 725 726 ! Exit loop and stop processing, if the last (requested) frame was encountered 727 IF (nframe >= last_frame) THEN 728 nframe = nframe - 1 729 CLOSE (UNIT=input_unit) 730 EXIT dcd_file_loop 731 END IF 732 733 END DO frame_loop 734 735 nframe = nframe - 1 736 737 CLOSE (UNIT=input_unit) 738 739 IF (iarg >= narg) EXIT dcd_file_loop 740 741 END DO dcd_file_loop 742 743 IF (info) THEN 744 WRITE (UNIT=error_unit,FMT="(A,/,A,I0)")& 745 "#",& 746 "# Frames processed: ",nframe_read 747 END IF 748 749 IF (ekin.AND.TRIM(ADJUSTL(id_dcd)) == "VEL") THEN 750 IF (info) THEN 751 WRITE (UNIT=error_unit,FMT="(A,/,A,F12.3)")& 752 "#",& 753 "# T [K] all frames: ",tavg/REAL(nframe_read,KIND=dp) 754 END IF 755 END IF 756 757 ! Print version information 758 IF (info) THEN 759 WRITE (UNIT=error_unit,FMT="(A)")& 760 "#",& 761 "# Normal termination of "//TRIM(version_info) 762 END IF 763 764 ! Close files 765 IF (have_atomic_labels) CLOSE (UNIT=xyz_input_unit) 766 IF (LEN_TRIM(out_file_name) > 0) CLOSE (UNIT=output_unit) 767 768 ! Cleanup 769 IF (ALLOCATED(atomic_label)) THEN 770 DEALLOCATE (atomic_label,STAT=istat) 771 IF (istat /= 0) THEN 772 CALL abort_program(routine_name,"Deallocation of the vector atomic_label failed") 773 END IF 774 END IF 775 776 IF (ALLOCATED(atomic_displacement)) THEN 777 DEALLOCATE (atomic_displacement,STAT=istat) 778 IF (istat /= 0) THEN 779 CALL abort_program(routine_name,"Deallocation of the vector atomic_displacement failed") 780 END IF 781 END IF 782 783 IF (ALLOCATED(atomic_mass)) THEN 784 DEALLOCATE (atomic_mass,STAT=istat) 785 IF (istat /= 0) THEN 786 CALL abort_program(routine_name,"Deallocation of the vector atomic_mass failed") 787 END IF 788 END IF 789 790 IF (ALLOCATED(atomic_temperature)) THEN 791 DEALLOCATE (atomic_temperature,STAT=istat) 792 IF (istat /= 0) THEN 793 CALL abort_program(routine_name,"Deallocation of the vector atomic_temperature failed") 794 END IF 795 END IF 796 797 IF (ALLOCATED(r)) THEN 798 DEALLOCATE (r,STAT=istat) 799 IF (istat /= 0) THEN 800 CALL abort_program(routine_name,"Deallocation of the array r failed") 801 END IF 802 END IF 803 804 IF (ALLOCATED(r0)) THEN 805 DEALLOCATE (r0,STAT=istat) 806 IF (istat /= 0) THEN 807 CALL abort_program(routine_name,"Deallocation of the array r0 failed") 808 END IF 809 END IF 810 811 IF (ALLOCATED(r_pbc)) THEN 812 DEALLOCATE (r_pbc,STAT=istat) 813 IF (istat /= 0) THEN 814 CALL abort_program(routine_name,"Deallocation of the array r_pbc failed") 815 END IF 816 END IF 817 818 IF (ALLOCATED(s)) THEN 819 DEALLOCATE (s,STAT=istat) 820 IF (istat /= 0) THEN 821 CALL abort_program(routine_name,"Deallocation of the array s failed") 822 END IF 823 END IF 824 825 IF (ALLOCATED(s_pbc)) THEN 826 DEALLOCATE (s_pbc,STAT=istat) 827 IF (istat /= 0) THEN 828 CALL abort_program(routine_name,"Deallocation of the array s_pbc failed") 829 END IF 830 END IF 831 832CONTAINS 833 834 SUBROUTINE abort_program(routine,message) 835 ! Abort the program after printing an error message to stderr 836 837 IMPLICIT NONE 838 839 CHARACTER(LEN=*), INTENT(IN) :: routine,message 840 841 CHARACTER(LEN=2*default_string_length) :: error_message 842 843 error_message = "*** ERROR in "//TRIM(routine)//": "//TRIM(message)//" ***" 844 WRITE (UNIT=default_error_unit,FMT="(/,A,/)") TRIM(error_message) 845 STOP "*** ABNORMAL PROGRAM TERMINATION of dumpdcd v3.0 ***" 846 847 END SUBROUTINE abort_program 848 849 SUBROUTINE build_h_matrix(a,b,c,alpha,beta,gamma,h) 850 ! Calculate the h matrix of a simulation cell given the cell edge lengths a, b, 851 ! and c in Angstrom and the angles alpha (<b,c)), beta (<(a,c)), and gamma (<(a,b)) 852 ! in degree. 853 854 IMPLICIT NONE 855 856 CHARACTER(LEN=*), PARAMETER :: routine_name = "build_h_matrix" 857 858 REAL(KIND=dp), INTENT(IN) :: a,alpha,b,beta,c,gamma 859 REAL(KIND=dp), DIMENSION(3,3), INTENT(OUT) :: h 860 861 REAL(KIND=dp) :: cosa,cosb,cosg,sing 862 863 cosa = COS(alpha/degree) 864 IF (ABS(cosa) < EPSILON(0.0_dp)) cosa = 0.0_dp 865 866 cosb = COS(beta/degree) 867 IF (ABS(cosb) < EPSILON(0.0_dp)) cosb = 0.0_dp 868 869 cosg = COS(gamma/degree) 870 IF (ABS(cosg) < EPSILON(0.0_dp)) cosg = 0.0_dp 871 872 sing = SIN(gamma/degree) 873 IF (ABS(sing) < EPSILON(0.0_dp)) sing = 0.0_dp 874 875 h(1,1) = 1.0_dp 876 h(2,1) = 0.0_dp 877 h(3,1) = 0.0_dp 878 879 h(1,2) = cosg 880 h(2,2) = sing 881 h(3,2) = 0.0_dp 882 883 h(1,3) = cosb 884 h(2,3) = (cosa - cosg*cosb)/sing 885 IF ((1.0_dp - h(1,3)**2 - h(2,3)**2) < 0.0_dp) THEN 886 CALL abort_program(routine_name,"Build of the h matrix failed, check cell information") 887 END IF 888 h(3,3) = SQRT(1.0_dp - h(1,3)**2 - h(2,3)**2) 889 890 h(:,1) = a*h(:,1) 891 h(:,2) = b*h(:,2) 892 h(:,3) = c*h(:,3) 893 894 END SUBROUTINE build_h_matrix 895 896 FUNCTION det_3x3(a) RESULT(det_a) 897 ! Returns the determinante of the 3x3 matrix a. 898 899 IMPLICIT NONE 900 901 REAL(KIND=dp), DIMENSION(3,3), INTENT(IN) :: a 902 903 REAL(KIND=dp) :: det_a 904 905 det_a = a(1,1)*(a(2,2)*a(3,3) - a(2,3)*a(3,2)) +& 906 a(1,2)*(a(2,3)*a(3,1) - a(2,1)*a(3,3)) +& 907 a(1,3)*(a(2,1)*a(3,2) - a(2,2)*a(3,1)) 908 909 END FUNCTION det_3x3 910 911 FUNCTION get_atomic_mass(element_symbol) RESULT(amass) 912 ! Get the atomic mass amass for an element 913 914 IMPLICIT NONE 915 916 CHARACTER(LEN=*), INTENT(IN) :: element_symbol 917 918 CHARACTER(LEN=*), PARAMETER :: routine_name = "get_atomic_mass" 919 920 REAL(KIND=dp) :: amass 921 922 SELECT CASE (TRIM(element_symbol)) 923 CASE ("O ") 924 amass = 15.9994_dp 925 CASE ("U ") 926 amass = 238.02891_dp 927 CASE DEFAULT 928 CALL abort_program(routine_name,"Unknown element symbol found") 929 END SELECT 930 931 amass = amass*massunit 932 933 END FUNCTION get_atomic_mass 934 935 SUBROUTINE invert_matrix_3x3(h,hinv,deth) 936 ! Calculate the inverse hinv and the determinant deth of the 3x3 matrix h. 937 938 IMPLICIT NONE 939 940 REAL(KIND=dp), DIMENSION(3,3), INTENT(IN) :: h 941 REAL(KIND=dp), DIMENSION(3,3), INTENT(OUT) :: hinv 942 REAL(KIND=dp), INTENT(OUT) :: deth 943 944 CHARACTER(LEN=*), PARAMETER :: routine_name = "invert_matrix_3x3" 945 946 deth = det_3x3(h) 947 948 ! Numerics 949 deth = ABS(deth) 950 IF (deth < 1.0E-10_dp) THEN 951 CALL abort_program(routine_name,"Invalid h matrix for cell found; det(h) < 1.0E-10") 952 END IF 953 954 hinv(1,1) = (h(2,2)*h(3,3) - h(3,2)*h(2,3))/deth 955 hinv(2,1) = (h(2,3)*h(3,1) - h(3,3)*h(2,1))/deth 956 hinv(3,1) = (h(2,1)*h(3,2) - h(3,1)*h(2,2))/deth 957 958 hinv(1,2) = (h(1,3)*h(3,2) - h(3,3)*h(1,2))/deth 959 hinv(2,2) = (h(1,1)*h(3,3) - h(3,1)*h(1,3))/deth 960 hinv(3,2) = (h(1,2)*h(3,1) - h(3,2)*h(1,1))/deth 961 962 hinv(1,3) = (h(1,2)*h(2,3) - h(2,2)*h(1,3))/deth 963 hinv(2,3) = (h(1,3)*h(2,1) - h(2,3)*h(1,1))/deth 964 hinv(3,3) = (h(1,1)*h(2,2) - h(2,1)*h(1,2))/deth 965 966 END SUBROUTINE invert_matrix_3x3 967 968 SUBROUTINE lowercase(string) 969 ! Convert all letters in a string to lowercase 970 CHARACTER(LEN=*), INTENT(INOUT) :: string 971 972 INTEGER :: i, iascii 973 974 DO i=1,LEN_TRIM(string) 975 iascii = ICHAR(string(i:i)) 976 IF ((iascii >= 65).AND.(iascii <= 90)) THEN 977 string(i:i) = CHAR(iascii + 32) 978 END IF 979 END DO 980 981 END SUBROUTINE lowercase 982 983 SUBROUTINE pbc(r,r_pbc,s,s_pbc,a,b,c,alpha,beta,gamma,debug,info,pbc0,h,hinv) 984 ! Apply the periodic boundary conditions (PBC) to the Cartesian coordinate array 985 ! r given the cell edge lengths a, b, and c in Angstrom and the angles alpha (<b,c)), 986 ! beta (<(a,c)), and gamma (<(a,b)) in degree. 987 ! On output r_pbc is updated with the "PBCed" input coordinates, s with the scaled 988 ! input coordinates, and s_pbc with scaled "PBCed" coordinates. 989 ! If pbc0 is true then fold to the range [-l/2,+l/2[ (origin at box centre) else fold 990 ! to the range [0,l[ (origin at lower left box corner). 991 992 IMPLICIT NONE 993 994 REAL(KIND=sp), DIMENSION(:,:), INTENT(IN) :: r 995 REAL(KIND=dp), DIMENSION(:,:), INTENT(OUT) :: r_pbc,s,s_pbc 996 REAL(KIND=dp), INTENT(IN) :: a,alpha,b,beta,c,gamma 997 LOGICAL, INTENT(IN) :: debug,info,pbc0 998 REAL(KIND=dp), DIMENSION(3,3), INTENT(OUT) :: h,hinv 999 1000 CHARACTER(LEN=*), PARAMETER :: routine_name = "pbc" 1001 1002 REAL(KIND=dp) :: deth 1003 INTEGER :: i,natom 1004 LOGICAL :: orthorhombic 1005 1006 natom= SIZE(r,1) 1007 IF (SIZE(r,2) /= 3) CALL abort_program(routine_name,"Array dimension for r must be 3") 1008 1009 ! Build h matrix 1010 1011 h(:,:) = 0.0_dp 1012 1013 IF ((alpha == 90.0_dp).AND.(beta == 90.0_dp).AND.(gamma == 90.0_dp)) THEN 1014 orthorhombic = .TRUE. 1015 h(1,1) = a 1016 h(2,2) = b 1017 h(3,3) = c 1018 ELSE 1019 orthorhombic = .FALSE. 1020 CALL build_h_matrix(a,b,c,alpha,beta,gamma,h) 1021 END IF 1022 1023 ! Build inverse of h 1024 hinv(:,:) = 0.0_dp 1025 CALL invert_matrix_3x3(h,hinv,deth) 1026 1027 IF (info) THEN 1028 WRITE (UNIT=error_unit,FMT="(A)") "#" 1029 IF (orthorhombic) THEN 1030 WRITE (UNIT=error_unit,FMT="(A)") "# Cell symmetry : orthorhombic" 1031 ELSE 1032 WRITE (UNIT=error_unit,FMT="(A)") "# Cell symmetry : non-orthorhombic" 1033 END IF 1034 IF (debug) THEN 1035 WRITE (UNIT=error_unit,FMT="(A)") "#" 1036 WRITE (UNIT=error_unit,FMT="(A,3F12.6,A)") "# / ",h(1,:)," \" 1037 WRITE (UNIT=error_unit,FMT="(A,3F12.6,A)") "# h = | ",h(2,:)," |" 1038 WRITE (UNIT=error_unit,FMT="(A,3F12.6,A)") "# \ ",h(3,:)," /" 1039 WRITE (UNIT=error_unit,FMT="(A)") "#" 1040 WRITE (UNIT=error_unit,FMT="(A,3F12.6,A)") "# / ",hinv(1,:)," \" 1041 WRITE (UNIT=error_unit,FMT="(A,3F12.6,A)") "# Inv(h) = | ",hinv(2,:)," |" 1042 WRITE (UNIT=error_unit,FMT="(A,3F12.6,A)") "# \ ",hinv(3,:)," /" 1043 WRITE (UNIT=error_unit,FMT="(A)") "#" 1044 WRITE (UNIT=error_unit,FMT="(A,F0.6)") "# det(h) = ",deth 1045 END IF 1046 END IF 1047 1048 ! Calculate scaled coordinates and wrap back all atoms, i.e. apply the PBC 1049 IF (orthorhombic) THEN 1050 ! Try to save some flops in the case of an orthorhombic box 1051 DO i=1,3 1052 s(:,i) = r(:,i)*hinv(i,i) 1053 END DO 1054 ELSE 1055 s(:,:) = MATMUL(r(:,:),TRANSPOSE(hinv(:,:))) 1056 END IF 1057 1058 IF (pbc0) THEN 1059 s_pbc(:,:) = s(:,:) - ANINT(s(:,:)) 1060 ELSE 1061 s_pbc(:,:) = s(:,:) - FLOOR(s(:,:)) 1062 END IF 1063 1064 IF (orthorhombic) THEN 1065 DO i=1,3 1066 r_pbc(:,i) = s_pbc(:,i)*h(i,i) 1067 END DO 1068 ELSE 1069 r_pbc(:,:) = MATMUL(s_pbc(:,:),TRANSPOSE(h(:,:))) 1070 END IF 1071 1072 END SUBROUTINE pbc 1073 1074 SUBROUTINE print_help() 1075 ! Print the program flags for help 1076 1077 IMPLICIT NONE 1078 1079 WRITE (UNIT=*,FMT="(T2,A)")& 1080 "",& 1081 "Program flags for "//TRIM(version_info)//":",& 1082 "",& 1083 " -debug, -d : Print debug information",& 1084 " -ekin : Dump just the ""temperature"" of each atom",& 1085 " -eformat : Print coordinates in scientific format",& 1086 " -eo : Write standard output and standard error to the same logical unit",& 1087 " -first_frame, -ff <int> : Number of the first frame which is dumped",& 1088 " -help, -h : Print this information",& 1089 " -info, -i : Print additional information for each frame (see also -debug flag)",& 1090 " -last_frame, -lf <int> : Number of the last frame which is dumped",& 1091 " -output, -o <file_name> : Name of the output file (default is stdout)",& 1092 " -output_format, -of <DCD|XMOL> : Output format for dump",& 1093 " -pbc : Apply the periodic boundary conditions (PBC) to each frame before it is dumped",& 1094 " (origin at lower left)",& 1095 " -pbc0 : Apply the periodic boundary conditions (PBC) to each frame before it is dumped",& 1096 " (origin at box centre)",& 1097 " -stride <int> : Stride for frame dump (allows to skip frames, e.g. by dumping each 10th frame)",& 1098 " -trace_atoms <real> : Print the atoms which left the simulation box given a threshold value in scaled units",& 1099 " -vel2cord, -v2c : Dump a VELocity DCD file as COoRDinate DCD file (hack which allows a digest by VMD)",& 1100 " -xyz_file, -xyz <file_name> : Name of a reference XYZ file in XMOL format that provides the atomic labels",& 1101 "" 1102 1103 WRITE (UNIT=*,FMT="(T2,A)")& 1104 "Usage examples:",& 1105 "",& 1106 " dumpdcd <optional flags> <DCD file(s)>",& 1107 "",& 1108 "Specific usage examples:",& 1109 "",& 1110 " dumpdcd UO2-2x2x2-pos-1.dcd (without atomic labels from XYZ file)",& 1111 " dumpdcd -xyz UO2-2x2x2.xyz UO2-2x2x2-pos-1.dcd (single DCD file)",& 1112 " dumpdcd -xyz UO2-2x2x2.xyz UO2-2x2x2-pos-1.dcd UO2-2x2x2-pos-2.dcd ... (multiple DCD files are dumped consecutively)",& 1113 " dumpdcd -info -xyz UO2-2x2x2.xyz UO2-2x2x2-pos-1.dcd UO2-2x2x2-pos-2.dcd (print additional information)",& 1114 " dumpdcd -debug -xyz UO2-2x2x2.xyz UO2-2x2x2-pos-1.dcd UO2-2x2x2-pos-2.dcd (print debug information)",& 1115 " dumpdcd -ekin -d -xyz UO2-2x2x2.xyz UO2-2x2x2-vel-1.dcd (print the ""temperature"" of each atom)",& 1116 " dumpdcd -ekin -xyz UO2-2x2x2.xyz UO2-2x2x2-vel-1.dcd (print just the temperature of each atom)",& 1117 " dumpdcd -first_frame 5 -last_frame 10 UO2-2x2x2-pos-1.dcd (just dump frame 5 to 10, ie. 6 frames in total)",& 1118 " dumpdcd -o outfile.xyz UO2-2x2x2-pos-1.dcd (write output to the file ""outfile.xyz"" instead of stdout)",& 1119 " dumpdcd -o test.xyz -output_format xmol -xyz ref.xyz -first 10 -last 10 test.dcd (dump 10th frame in XMOL format)",& 1120 " 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)",& 1121 " 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)",& 1122 " dumpdcd -o part.dcd -of dcd -first 10 -lf 100 -stride 10 test.dcd (dump the frames 10, 20, ..., 100 to the file part.dcd)",& 1123 " dumpdcd -output new.dcd -output_format dcd -pbc old.dcd (dump all frames applying PBC to the output file new.dcd)",& 1124 " dumpdcd -o new.dcd -of dcd -pbc -trace_atoms 0.02 old.dcd (all atoms more than 2% out of the box are listed)",& 1125 " dumpdcd -o new.dcd -e out_of_box.log -of dcd -pbc -trace_atoms 0.1 old.dcd (atoms more than 10% out of the box are listed)",& 1126 " dumpdcd -o new.dcd -of dcd -vel2cord old.dcd (dump old.dcd as new.dcd and change only the DCD id string from VEL to CORD)",& 1127 " dumpdcd -o new.dcd -of dcd -frc2cord old.dcd (dump old.dcd as new.dcd and change only the DCD id string from FRC to CORD)",& 1128 " dumpdcd -i -disp UO2-2x2x2-pos-1.dcd (dump the displacements of all atoms w.r.t. their positions in the first frame)",& 1129 " dumpdcd -i -of dcd -disp UO2-2x2x2-pos-1.dcd (dump the atomic displacements as x-coordinates of a DCD CORD file)",& 1130 " dumpdcd -i -of dcd -ekin -v2c -xyz UO2-2x2x2.xyz UO2-2x2x2-vel-1.dcd (dump the atomic temperatures as x-coordinates of a ",& 1131 " DCD CORD file -> hack for VMD)",& 1132 "" 1133 1134 WRITE (UNIT=*,FMT="(T2,A)")& 1135 "Notes:",& 1136 "",& 1137 " - For -ekin a XYZ file is required to obtain the atomic labels",& 1138 " - The -info and the -debug flags provide a more detailed output which is especially handy for tracing problems",& 1139 " - The output in DCD format is in binary format",& 1140 " - The input coordinates should be in Angstrom. Velocities and forces are expected to be in atomic units",& 1141 "" 1142 1143 END SUBROUTINE print_help 1144 1145 SUBROUTINE uppercase(string) 1146 ! Convert all letters in a string to uppercase 1147 1148 IMPLICIT NONE 1149 1150 CHARACTER(LEN=*), INTENT(INOUT) :: string 1151 1152 INTEGER :: i,iascii 1153 1154 DO i=1,LEN_TRIM(string) 1155 iascii = ICHAR(string(i:i)) 1156 IF ((iascii >= 97).AND.(iascii <= 122)) THEN 1157 string(i:i) = CHAR(iascii - 32) 1158 END IF 1159 END DO 1160 1161 END SUBROUTINE uppercase 1162 1163 SUBROUTINE write_out_of_box_atoms(atomic_label,r,s,eps_out_of_box,h) 1164 ! Print a list of all atoms which have left the simulation box 1165 1166 IMPLICIT NONE 1167 1168 CHARACTER(LEN=5), DIMENSION(:), INTENT(IN) :: atomic_label 1169 REAL(KIND=sp), DIMENSION(:,:), INTENT(IN) :: r 1170 REAL(KIND=dp), DIMENSION(:,:), INTENT(IN) :: s 1171 REAL(KIND=dp), INTENT(IN) :: eps_out_of_box 1172 REAL(KIND=dp), DIMENSION(3,3), INTENT(IN) :: h 1173 1174 REAL(KIND=dp), DIMENSION(3) :: ds,dr 1175 REAL(KIND=dp) :: rl,s_max,s_min,sl 1176 INTEGER :: i,iatom,natom,ncount 1177 1178 ! Quick return, if no action is requested 1179 IF (eps_out_of_box <= 0.0_dp) RETURN 1180 1181 s_max = 1.0_dp + eps_out_of_box 1182 s_min = -eps_out_of_box 1183 natom = SIZE(s,1) 1184 ncount = 0 1185 DO iatom=1,natom 1186 IF (ANY(s(iatom,:) < s_min).OR.& 1187 ANY(s(iatom,:) > s_max)) THEN 1188 ncount = ncount + 1 1189 IF (ncount == 1) THEN 1190 WRITE (UNIT=error_unit,FMT="(A)")& 1191 "#",& 1192 "# Atoms out of box:",& 1193 "# Atom index label x y z |dr| |ds|" 1194 END IF 1195 ds(:) = s(iatom,:) 1196 DO i=1,3 1197 IF (s(iatom,i) < 0.0_dp) ds(i) = 0.0_dp 1198 IF (s(iatom,i) >= 1.0_dp) ds(i) = 1.0_dp 1199 END DO 1200 ds(:) = s(iatom,:) - ds(:) 1201 sl = SQRT(ds(1)**2 + ds(2)**2 + ds(3)**2) 1202 dr(:) = MATMUL(h(:,:),ds(:)) 1203 rl = SQRT(dr(1)**2 + dr(2)**2 + dr(3)**2) 1204 WRITE (UNIT=error_unit,FMT="(A,I10,1X,A5,5(1X,F14.6))")& 1205 "# ",iatom,ADJUSTR(atomic_label(iatom)),r(iatom,:),rl,sl 1206 END IF 1207 END DO 1208 WRITE (UNIT=error_unit,FMT="(A,I0,A)") "# ",ncount," atom(s) out of box" 1209 1210 END SUBROUTINE write_out_of_box_atoms 1211 1212END PROGRAM dumpdcd 1213