1! ------------------------------------------------------------- 2! Copyright (c) Universite de Montreal 3! M-A.Malouin and F.El-Mellouhi, 2005 4! 5! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 6! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 7! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 8! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 9! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 10! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 11! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 12! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 13! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 14! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 15! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 16! 17! ------------------------------------------------------------- 18! Program use to write OpenDx file in .dx format from files 19! produce by siesta runs (use in VEP OpenDx program 'Dx View.net') 20! 21! Write data for: 1) Electronic density 22! 2) Atoms 23! 3) Lattice Vectors 24! 4) Unit Cell Frame 25! 26! Written in fortran 90/95 ... use the same compiler as with siesta 27! for correct reading of unformatted binary file 28! 29! The following compilation has been working fine...on linux and mac_os_x 30! e.g: xlf95 (or xlf90) -o DxFormat DxFormat.f 31! ifort -o DxFormat DxFormat.f90 (you must change the file extension for this one) 32! 33! ------------------------------------------------------------- 34! This program convert siesta output from 35! .RHO, .DRHO, .VH, .VT, .IOCH, .TOCH, .LDOS and .XV files 36! in .dx format and combine them in one big file for use 37! with OpenDx visualization program DxView.net 38! 39! Notes on options: 40! 1) Siesta uses periodic cell representation and some 41! shifting corrections must in general be done for 42! isosurfaces to ensure centered visualization 43! 2) Atoms and isosurface are not writen in the same 44! coordinates system by Siesta and translation must 45! be done to ensure consistency with isosurface 46! (Unfortunately this work only for orthogonal cell 47! right now...but this will eventually be fixed) 48! 3) When comparing your results with your initial 49! configuration, viewing only change data can be 50! useful, so you can substract from your data any 51! configuration that you want to use as a base 52! reference to view only the changes 53! (Note: The difference is applied before conversion 54! to .dx file format for OpenDx) 55! 56! Special note: If the file systemLabel.xyz is found, the 57! legend labels will be replaced by their 58! real labels instead of their default id code 59! 60! Usage: Dxformat system_label [-s, -t, -d diff_label] 61! system_label: the one defined in the .fdf input file used 62! by Siesta for your simulation 63! Options: -s --> shifting of coordinates for isosurfaces (1) 64! -t --> translation of coordinates for atoms (2) 65! -d --> use files from second system_label 'diff_label' 66! as Reference base data (3) 67! 68! PLEASE READ THE README FILE FOR MORE INFORMATIONS ABOUT USING THIS 69! PROGRAM TOGETHER WITH SIESTA AND OPENDX 70! 71! Version 1.0 72! @author Marc-Andre Malouin 73! @date July 12, 2005 74! ------------------------------------------------------------- 75 Program DxFormat 76 Implicit none 77 78 integer, parameter :: sp = selected_real_kind(6,30) 79 integer, parameter :: dp = selected_real_kind(14,100) 80 81 Integer, Parameter :: MAX_LENGTH = 80, EXT_LENGTH = 5, OP_LENGTH = 2 ! For character array variables 82 Integer, Parameter :: READ_UNIT = 1, WRITE_UNIT = 2, DIFF_UNIT = 3, XYZ_UNIT = 4 ! IO/Units 83 84! Index for loops 85 Integer :: i 86 87! To check for i/o and allocation error 88 Integer :: status = 0 ! If not zero after use: an error occured ! 89 90! Various file names and extensions 91 Character(MAX_LENGTH) :: system_label, diff_label, filename 92! The different type of isosurface process by this program 93 Character(EXT_LENGTH), Parameter :: FILE_FORMAT(8)=(/".RHO ",".DRHO", ".VH ", ".VT ",".IOCH",".TOCH",".LDOS",".XV "/) 94 Character(EXT_LENGTH), Parameter :: ISO = ".ISO ", DX = ".dx ", XYZ = ".xyz " 95 96 Character(6), Parameter :: NO_ISO = '"None"' 97 Integer :: nb_iso = 0, xv = 0 98 Logical :: is_fileFound(7)= (/.false.,.false.,.false.,.false., .false.,.false.,.false./) 99 100! To write Unit Cell and Lattice Vectors info only once for iso files 101 Logical :: add_uclv = .true. 102 103! Command line option input buffer 104 Character(MAX_LENGTH) :: buffer 105 106! Command line arguments 107 Integer :: nargs ! Number of arguments 108 Integer, External :: iargc ! Function to retreive nargs 109 110! Valid command line options 111 Character(OP_LENGTH), Parameter :: OPT_SHIFT_ISO = "-s" 112 Character(OP_LENGTH), Parameter :: OPT_TRANSLATE_ATOMS = "-t" 113 Character(OP_LENGTH), Parameter :: OPT_DIFF_ISO = "-d" 114 115! Default options behavior 116 Logical :: shift_iso = .false. 117 Logical :: translate_atoms = .false. 118 Logical :: diff = .false. 119 120! --- Io/Formats--- 121 10 Format ("Option '",A, "' use : ",A) ! For translation option 122 20 Format ('Bad option argument: ',A) ! For bad argument 123 124! ------------------------------------------------------------- 125! ---Reading command line arguments--- 126! ------------------------------------------------------------- 127 128! Number of arguments 129 nargs = iargc() 130 If(nargs == 0) Then 131 Call printInfo() 132 Call printUsage() 133 Stop 134 End if 135 136! First argument require: the system_label 137 Call getarg(1, system_label) 138 Write(*,*) "Processing for system_label: ", system_label 139 140! No difference by default (difference only if diff_label != system_label) 141 diff_label = system_label 142 143! Reading optional arguments 144 i=2 145 Do while (i <= nargs) 146 Call getarg(i, buffer) 147 148 Select case (buffer) 149 Case (OPT_SHIFT_ISO) 150 Write(*,10) trim(buffer), 'shifting of coordinates for isosurfaces' 151 shift_iso = .true. 152 Case (OPT_TRANSLATE_ATOMS) 153 Write(*,10) trim(buffer), 'translation of coordinates for atoms' 154 translate_atoms = .true. 155 Case (OPT_DIFF_ISO) 156 i=i+1 157 Call getarg(i, diff_label) 158 159 If(trim(diff_label) == trim(system_label)) Then 160 Write(*,10) trim(buffer), 'BUT diff_label = system_label...option skip (not used)' 161 Else 162 Write(*,10) trim(buffer), "using files from system_label '" // trim(diff_label) // "' as reference for difference" 163 End if 164 Case DEFAULT 165 Write(*,20) trim(buffer) 166 Call printUsage() 167 Stop 168 End select 169 i=i+1 170 End do 171 172! ------------------------------------------------------------- 173! ---Processing files--- 174! ------------------------------------------------------------- 175 176 filename = trim(system_label)//trim(DX) 177 Open(UNIT=WRITE_UNIT, FILE=trim(filename), STATUS='Unknown', ACTION='write', IOSTAT=status, FORM='formatted') 178 179! IOError checking 180 If(status /= 0) Then 181 Write(*,*) 'Error while creating or replacing the file: ' // trim(filename) 182 Stop ! End of program 183 End if 184 185 Do i=1, size(FILE_FORMAT) 186 187! Opening file 188 filename = trim(system_label)//trim(FILE_FORMAT(i)) 189 Write(*,*) 190 Write(*,*) 'Reading data from file: ', trim(filename) 191 192 If(i == size(FILE_FORMAT)) Then 193 Open(UNIT=READ_UNIT,FILE=filename,STATUS='old',ACTION='READ',IOSTAT=status,FORM='formatted') 194 Else 195 Open(UNIT=READ_UNIT,FILE=filename,STATUS='old',ACTION='READ', IOSTAT=status,FORM='unformatted') 196 End if 197 198 If(status == 0) Then 199 is_fileFound(i) = .true. 200 201 If(trim(diff_label) /= trim(system_label)) Then 202 diff = .true. 203! Opening file for difference option 204 filename = trim(diff_label)//trim(FILE_FORMAT(i)) 205 206 If(i == size(FILE_FORMAT)) Then 207 Open(UNIT=DIFF_UNIT,FILE=filename,STATUS='old',ACTION='READ',IOSTAT=status,FORM='formatted') 208 Else 209 Open(UNIT=DIFF_UNIT,FILE=filename,STATUS='old',ACTION='READ',IOSTAT=status,FORM='unformatted') 210 End if 211 212 If(status /= 0) Then 213 Write(*,*) 'File "' // trim(filename) // '" for difference not found...option not used' 214 diff = .false. ! Option not used 215 End if 216 End if 217 218! Processing files 219 If(i == size(FILE_FORMAT)) Then 220 221! Is the systemLabel.xyz present ? 222 filename = trim(system_label)//trim(XYZ) 223 Open(UNIT=XYZ_UNIT,FILE=filename,STATUS='old',ACTION='READ',IOSTAT=status,FORM='formatted') 224 225 If(status == 0) Then 226 Write(*,*) '"' // trim(filename) // '" file found: changing legend labels to strings instead of numbers' 227 Call processAtoms(FILE_FORMAT(i), translate_atoms, diff,.true.) ! Yes 228 Else 229 Call processAtoms(FILE_FORMAT(i), translate_atoms, diff,.false.) ! No 230 End if 231 232 xv = 1 233 Close(XYZ_UNIT) 234 235 Else 236 Call processIso(FILE_FORMAT(i), shift_iso, diff, add_uclv) 237 add_uclv = .false. 238 nb_iso = nb_iso + 1 239 End if 240 241 If(diff) Close(DIFF_UNIT) 242 243 Else ! An error occured or the file was not found 244 Write(*,*) 'Error while opening the file: ' // trim(filename) // ' --> Skipping file...' 245 End if 246 247 Close(READ_UNIT) 248 249 End do 250 251! Grouping members and data together in .dx file 252 Write(WRITE_UNIT,*) '############################################' 253 Write(WRITE_UNIT,*) '#' 254 Write(WRITE_UNIT,*) '# GROUPING MEMBERS ' 255 Write(WRITE_UNIT,*) '#' 256 Write(WRITE_UNIT,*) '############################################' 257 Write(WRITE_UNIT,*) 'object "iso" array type string rank 1 shape ',(EXT_LENGTH+1),' items ', nb_iso+1, ' data follows' 258 Write(WRITE_UNIT,*) NO_ISO 259 Do i=1, size(is_fileFound) 260 If(is_fileFound(i)) Write(WRITE_UNIT,*) '"',trim(FILE_FORMAT(i)),'"' 261 End do 262 Write(WRITE_UNIT,*) 'attribute "XV" value number ', xv 263 Write(WRITE_UNIT,*) '############################################' 264 Write(WRITE_UNIT,*) 'object "all" class group' 265 Write(WRITE_UNIT,*) 'member 0 value "iso"' 266 Do i=1, size(is_fileFound) 267 If(is_fileFound(i)) Write(WRITE_UNIT,*) 'member "' // trim(FILE_FORMAT(i)) // '" value "iso' // trim(FILE_FORMAT(i)) // '"' 268 End do 269 If(xv == 1) Write(WRITE_UNIT,*) 'member "' // trim(FILE_FORMAT(size(FILE_FORMAT))) // '" value "molecule"' 270 If(nb_iso /= 0) Then 271 Write(WRITE_UNIT,*) 'member "vectors" value "vectors' // trim(ISO) // '"' 272 Write(WRITE_UNIT,*) 'member "cell" value "cell' // trim(ISO) // '"' 273 End if 274 275 close(WRITE_UNIT) 276 277 End program 278 279! ------------------------------------------------------------------ 280! Subroutines for processing files 281! ------------------------------------------------------------------ 282 Subroutine processIso (type, translate, use_diff, add_uclv) 283 Implicit none 284 285 integer, parameter :: sp = selected_real_kind(6,30) 286 integer, parameter :: dp = selected_real_kind(14,100) 287 288 Integer, Parameter :: EXT_LENGTH = 5 ! For type length 289 Integer, Parameter :: READ_UNIT = 1, WRITE_UNIT = 2, DIFF_UNIT = 3 ! IO/Units 290 Character(EXT_LENGTH), Parameter :: ISO = ".ISO " 291 292! Parameters 293 Character(EXT_LENGTH), Intent(IN) :: type ! type = isosurface type (RHO, DRHO, VH, VT, IOCH or TOCH) 294 Logical, Intent(IN) :: translate ! Translation correction ? 295 Logical, Intent(IN) :: add_uclv ! Write Cell And Vector info ? 296 Logical, Intent(IN) :: use_diff ! Use difference option or not 297 298 Logical :: diff 299 300! For storing isosurfaces data read from input_file 301! (N.B: We don't use nsd...just use to read an extra token in file...number spin density information) 302 303 Integer :: nsd, ntp ! ntp = ntm(1) * ntm(2) * ntm(3) (mesh volume) 304 Integer :: ntm(3), tntm(3) ! The mesh dimension in x, y, z 305 ! (tntm to test system compatibility when use_diff is used) 306 real(dp) :: ucell(3,3), tcell(3,3) ! Lattice unit cell 307 ! (tcell to test system compatibility when use_diff is used) 308 309 310! 3D Electronic density projected in one dimension (at each x,y,z of mesh) 311! [dt is after translation (shift of origin for correct visualization)] 312 real(sp), Allocatable :: ed(:), edt(:), ted(:) ! ted for difference option (when use_diff is specified) 313 314 Integer :: delta(3) ! To shift the origin by an amount delta in x, y, z 315 316! For allocation errors testing 317 Integer :: status 318 319! Dummy variable 320 Integer :: temp 321 322! Index for various loops (ix,iy,iz,iix,iiy,iiz are use for triple loop for translation) 323 Integer :: i, j, ix, iy, iz, iix, iiy, iiz 324 325 326! --- Io/Formats--- 327 200 Format (3F30.25) ! For printing matrix rows (lattice vectors) of ucell 328 210 Format (A, ' ---> ',3F30.25) 329 220 Format (3F15.10) ! A (x,y,z) point 330 331! --Reading-- 332 Read(READ_UNIT) ucell 333 Read(READ_UNIT) ntm, nsd 334 335 ntp = ntm(1) * ntm(2) * ntm(3) 336 337! Difference (Testing system compatibility) 338 If(use_diff) Then 339 diff = .true. 340 341 Read(DIFF_UNIT) tcell 342 Read(DIFF_UNIT) tntm, temp 343 344 If(tntm(1) /= ntm(1) .OR. tntm(2) /= ntm(2) .OR. tntm(3) /= ntm(3) .OR. temp /= nsd) Then 345 diff = .false. 346 Else 347 Do i=1, 3 348 If(tcell(i,1) /= ucell(i,1) .OR. tcell(i,2) /= ucell(i,2) .OR. tcell(i,3) /= ucell(i,3)) Then 349 diff = .false. 350 Exit 351 End if 352 End do 353 End if 354 355 If(.NOT. diff) Write(*,*) 'Incompatible system specified for difference...option not used' 356 Else 357 diff = .false. 358 End if 359 360! Printing values in terminal 361 If(add_uclv) Then 362 Write(*,*) 363 Write(*,*) 'Common data for all isosurface files:' 364 Write(*,*) 'Ntm size = ', ntp, ' = ', ntm(1), ' * ', ntm(2), ' * ', ntm(3), ' ( nsd = ', nsd,')' 365 Write(*,*) 'Cell Vectors:' 366 Write(*,210) "x'",(ucell(1,j), j=1, 3) 367 Write(*,210) "y'",(ucell(2,j), j=1, 3) 368 Write(*,210) "z'",(ucell(3,j), j=1, 3) 369 End if 370 371! Array allocation 372 status = 0 373 Allocate(ed(ntp),edt(ntp), STAT=status) 374 If(diff) Allocate(ted(ntm(1)), STAT=status) 375 If(status /= 0) Then ! Allocation error checking 376 Write(*,*) 'Array allocation error (overflow ???)...ending program !' 377 Stop 378 End if 379 380 Write(*,*) 381 Write(*,*) 'Reading charge density ... ' 382 Do iz=0, ntm(3)-1 383 Write(*,*) 'Reading plane ', iz+1, ' ...' 384 Do iy=0, ntm(2)-1 385 Read(READ_UNIT) (ed(ix+iy*ntm(1)+iz*ntm(1)*ntm(2)), ix=1,ntm(1)) 386! Difference 387 If(diff) Then 388 Read(DIFF_UNIT) (ted(i),i=1,ntm(1)) 389 Do i = 1, ntm(1) 390 ed(i+iy*ntm(1)+iz*ntm(1)*ntm(2)) = ed(i+iy*ntm(1)+iz*ntm(1)*ntm(2)) - ted(i) 391 End do 392 End if 393 End do 394 End do 395 396! Shifting of the origin for centered visualization of isosurface 397 If(translate) Then 398 399 Write(*,*) 400 Write(*,*) 'Shifting of coordinates for visualization ...' 401 Do i=1, 3 402 delta(i) = ntm(i)/2.0 ! ucell(i,i)/2.0 * ntm(i)/ucell(i,i) 403 End do 404 Write(*,*) 'Delta shifting coefficients in x, y, z:', delta 405 406 Do iz=1, ntm(3) 407 Do iy=1, ntm(2) 408 Do ix=1, ntm(1) 409 410! Wise translation for not going out of bounds 411 If(ix > ntm(1)/2) Then 412 iix = ix - delta(1) 413 Else if(ix < ntm(1)/2) Then 414 iix = ix + delta(1) 415 End if 416 If(iy > ntm(2)/2) Then 417 iiy = iy - delta(2) 418 Else if(iy < ntm(2)/2) Then 419 iiy = iy + delta(2) 420 End if 421 If(iz > ntm(3)/2) Then 422 iiz = iz - delta(3) 423 Else if(iz < ntm(3)/2) Then 424 iiz = iz + delta(3) 425 End if 426 427! Bounds testing (Just in case...should not happen !) 428 If(iix < 1) Stop 'iix < 0' 429 If(iiy < 1) Stop 'iiy < 0' 430 If(iiz < 1) Stop 'iiz < 0' 431 If(iix > ntm(1)) Stop 'iix > cell' 432 If(iiy > ntm(2)) Stop 'iiy > cell' 433 If(iiz > ntm(3)) Stop 'iiz > cell' 434 435! Translation of values around the new shifted origin 436 i = ix + (iy-1)*ntm(1) + (iz-1)*ntm(1)*ntm(2) 437 j = iix + (iiy-1)*ntm(1) + (iiz-1)*ntm(1)*ntm(2) 438 edt(j) = ed(i) 439 End do 440 End do 441 End do 442 443! Final step (copy of translated coordinates) 444 ed = edt 445 446 End if 447 448! Writing data in .dx format 449 450! Unit cell and lattice vectors info 451 If(add_uclv) Call write_UCLV(ucell, ISO) 452 453 Write(WRITE_UNIT,*) '############################################' 454 Write(WRITE_UNIT,*) '#' 455 Write(WRITE_UNIT,*) '# ' // trim(type) // 'FILE INFO' 456 Write(WRITE_UNIT,*) '#' 457 Write(WRITE_UNIT,*) '############################################' 458 459 Write(WRITE_UNIT,*) 'object "positions' // trim(type) // '" class gridpositions counts', ntm(3), ntm(2), ntm(1) 460 461 Write(WRITE_UNIT,*) 'origin 0 0 0' 462 Do i=3, 1, -1 463 Write(WRITE_UNIT,*) 'delta',((ucell(i,j)/ntm(i)),j=1,3) 464 End do 465 466 Write(WRITE_UNIT,*) 'object "connections' // trim(type) // '" class gridconnections counts', ntm(3), ntm(2), ntm(1) 467 Write(WRITE_UNIT,*) 'object "data' // trim(type) // '" class array type float rank 0 items', ntp, 'data follows' 468 469 Write(WRITE_UNIT,*) (ed(i), i=1, ntp) 470 471 Write(WRITE_UNIT,*) 'object "iso' // trim(type) // '" class field' 472 Write(WRITE_UNIT,*) 'component "positions" value "positions' // trim(type) // '"' 473 Write(WRITE_UNIT,*) 'component "connections" value "connections' // trim(type) // '"' 474 Write(WRITE_UNIT,*) 'component "data" value "data' // trim(type) // '"' 475 Write(WRITE_UNIT,*) 476 477 End subroutine 478! ------------------------------------------------------------------ 479 Subroutine processAtoms (type, translate, Use_diff, add_xyz) 480 Implicit none 481 482 integer, parameter :: sp = selected_real_kind(6,30) 483 integer, parameter :: dp = selected_real_kind(14,100) 484 485 Integer, Parameter :: MAX_LENGTH = 80, EXT_LENGTH = 5 ! For type length 486 Integer, Parameter :: READ_UNIT = 1, WRITE_UNIT = 2, DIFF_UNIT = 3, XYZ_UNIT = 4 ! IO/Units 487 488! Parameters 489 Character(EXT_LENGTH), Intent(IN) :: type ! type = atom type (.XV here) 490 Logical, Intent(IN) :: translate 491 Logical, Intent(IN) :: use_diff 492 Logical, Intent(IN) :: add_xyz 493 494 Logical :: diff, xyz 495 496! For storing atoms data from input file 497 real(dp) :: ucell(3,3), tcell(3,3) ! Lattice unit cell (tcell to test system compatibility when use_diff is used) 498 real(dp) :: N(3) !To translate the origin by an amount dx, dy and dz 499 Character(MAX_LENGTH), allocatable :: atom_label(:) ! Atom labels (code number or string if add_xyz specified) 500 Integer, allocatable :: atomic_number(:) ! Atoms atomic number 501 real(dp), allocatable :: x(:), y(:), z(:) ! Coordinates 502 real(dp) :: tx, ty, tz ! For difference opion (when use_diff is specified) 503 Integer :: nb_atoms ! Total number of atoms (future size of all allocatable variables above) 504 505 Integer :: ndiff_atoms ! Total number of different atoms 506 Logical, allocatable :: diff_atom(:) ! If the atom moved or not (Use only if use_diff is specified) 507 508! For allocation error testing 509 Integer :: status 510 511! Dummy variable 512 Integer :: temp 513 514! Index for various loops 515 Integer :: i, j 516 517! ---Io/Formats---- 518 200 Format (3F30.25) ! For printing matrix rows (lattice vectors) of ucell 519 210 Format (A, ' ---> ',3F30.25) 520 220 Format (3F15.10) ! A (x,y,z) point 521 522! --Reading-- 523 Do i=1, 3 524 Read(READ_UNIT,*) (ucell(i,j), j=1,3) 525 End do 526 Read(READ_UNIT,*) nb_atoms ! Number of atoms 527 528! Get .xyz file ready to read 529 If(add_xyz) Then 530 xyz = .true. 531 Read(XYZ_UNIT,*) temp ! Skipping first line (the numbers of atoms) 532 If(temp /= nb_atoms) Then 533 xyz = .false. 534 Write(*,*) 'Warning ! Number of atoms in file '//trim(type)//' and file .xyz do not match...legend labels not read...' 535 End if 536 Else 537 xyz = .false. 538 End if 539 540! Difference (Testing for system compatibility) 541 If(use_diff) Then 542 diff = .true. 543 Do i=1, 3 544 Read(DIFF_UNIT,*) (tcell(i,j), j=1,3) 545 End do 546 Read(DIFF_UNIT,*) temp ! Number of atoms 547 548 If(temp /= nb_atoms) Then 549 diff = .false. 550 Else 551 Do i=1, 3 552 If(tcell(i,1) /= ucell(i,1) .OR. tcell(i,2) /= ucell(i,2) .OR. tcell(i,3) /= ucell(i,3)) Then 553 diff = .false. 554 Exit 555 End if 556 End do 557 End If 558 559 If(.NOT. diff) Write(*,*) 'Incompatible system specified for difference...option not used' 560 Else 561 diff = .false. 562 End if 563 564! Printing of values in terminal 565 Write(*,*) 'Number of atoms = ', nb_atoms 566 Write(*,*) 'Cell Vectors:' 567 Write(*,210) "x'",(ucell(1,j), j=1, 3) 568 Write(*,210) "y'",(ucell(2,j), j=1, 3) 569 Write(*,210) "z'",(ucell(3,j), j=1, 3) 570 571! Array allocation 572 status = 0 573 Allocate(x(nb_atoms), y(nb_atoms), z(nb_atoms),STAT=status) 574 If(status /= 0) Then ! Allocation error checking 575 Write(*,*) 'Array allocation error (overflow ???)...ending program !' 576 Stop 577 End if 578 Allocate (atom_label(nb_atoms), atomic_number(nb_atoms), diff_atom(nb_atoms), STAT=status) 579 If(status /= 0) Then ! Allocation error checking 580 Write(*,*) 'Array allocation error (overflow ???)...ending program !' 581 Stop 582 End if 583 584 ndiff_atoms = 0 585 Do i=1,nb_atoms 586 587 Read(READ_UNIT,*)atom_label(i),atomic_number(i),x(i),y(i),z(i) 588 589! Replacing labels with text if .xyz file specified 590 If(xyz) Read(XYZ_UNIT,*) atom_label(i), tx, ty, tz ! Only first data will be use..skipping already read coordinates x, y, z 591 592 diff_atom(i) = .true. ! All atoms are different from nothing by default...unless otherwise specified 593 594! Difference 595 If(diff) Then 596 Read(DIFF_UNIT,*)temp,temp,tx,ty,tz 597! If the atom moved, we keep it (it is different)...otherwise we remove it 598 If(x(i) == tx .AND. y(i) == ty .AND. z(i) == tz) Then 599 diff_atom(i) = .false. 600 ndiff_atoms = ndiff_atoms + 1 601 End if 602 End if 603 End do 604 605! Last consistency test for difference option 606 If(nb_atoms == ndiff_atoms) Then 607 Write(*,*) 'File use for difference is identical to original file...option not used' 608 ndiff_atoms = 0 609 Do i=1, nb_atoms 610 diff_atom(i) = .true. 611 End do 612 End if 613 614! Translation of coordinates for consistency with isosurface 615 If(translate) Then 616 617 Do i=1, 3 618! Norm of lattice vectors 619 N(i) = sqrt(ucell(i,1)**2 + ucell(i,2)**2 + ucell(i,3)**2) 620 End do 621 622 Write(*,*) 623 Write(*,*) "Delta translation coefficients for atoms in each lattice vectors direction x', y', z': " 624 Write(*,200) N/2 625 626! --Translating atom coordinates-- 627 Do i=1, nb_atoms 628 x(i) = x(i) + (ucell(1,1) + ucell(2,1) + ucell(3,1))/2 629 y(i) = y(i) + (ucell(1,2) + ucell(2,2) + ucell(3,2))/2 630 z(i) = z(i) + (ucell(1,3) + ucell(2,3) + ucell(3,3))/2 631 End do 632 End if 633 634 Write(*,*) 635 Write(*,*) 'Writing atoms data...' 636 637! Writing to file 638 Write(WRITE_UNIT,*) '############################################' 639 Write(WRITE_UNIT,*) '#' 640 Write(WRITE_UNIT,*) '# ATOMS COORDINATES AND COLOR INFO' 641 Write(WRITE_UNIT,*) '#' 642 Write(WRITE_UNIT,*) '############################################' 643 644! Atoms positions 645 Write(WRITE_UNIT,*) 'object "atomcoord" array type float rank 1 shape 3 items ', (nb_atoms-ndiff_atoms), ' data follows' 646 Do i=1, nb_atoms 647 If(diff_atom(i)) Write(WRITE_UNIT,220) x(i), y(i), z(i) 648 End do 649 650! Atoms size 651 Write(*,*) ' Default atoms size: Base on unique atomic number...' 652 653 Write(WRITE_UNIT,*) 'object "atomsize" array type float rank 0 items ', (nb_atoms-ndiff_atoms), ' data follows' 654 Do i=1, nb_atoms 655 If(diff_atom(i)) Write(WRITE_UNIT,*) (atomic_number(i) + 0.0) 656 End do 657! Atoms code number id 658 Write(*,*) ' Default atoms color: Base on unique code number...' 659 660 Write(WRITE_UNIT,*) & 661 'object "label" array type string rank 1 shape ', & 662 MAX_LENGTH,' items ', (nb_atoms-ndiff_atoms), ' data follows' 663 Do i=1, nb_atoms 664 If(diff_atom(i)) Write(WRITE_UNIT,*) ('"'//trim(atom_label(i))//'"') 665 End do 666 Write(WRITE_UNIT,*) 'attribute "dep" value string "positions"' 667 668 Write(WRITE_UNIT,*) 'object "atoms' // trim(type) // '" class field' 669 Write(WRITE_UNIT,*) 'component "positions" value "atomcoord"' 670 Write(WRITE_UNIT,*) 'component "data" value "atomsize"' 671 Write(WRITE_UNIT,*) 'component "id" value "label"' 672 Write(WRITE_UNIT,*) 673 674! Unit cell and lattice vectors info 675 Call write_UCLV(ucell, type) 676 677 Write(WRITE_UNIT,*) '############################################' 678 Write(WRITE_UNIT,*) 'object "molecule" class group' 679 Write(WRITE_UNIT,*) 'member "atoms" value "atoms' // trim(type) // '"' 680 Write(WRITE_UNIT,*) 'member "vectors" value "vectors' // trim(type) // '"' 681 Write(WRITE_UNIT,*) 'member "cell" value "cell' // trim(type) // '"' 682 Write(WRITE_UNIT,*) 683 684 End subroutine 685! ------------------------------------------------------------------ 686 Subroutine write_UCLV (ucell, label) 687 Implicit none 688 689 integer, parameter :: sp = selected_real_kind(6,30) 690 integer, parameter :: dp = selected_real_kind(14,100) 691 692 Integer, Parameter :: WRITE_UNIT = 2 ! IO/Units 693 694! Parameters 695 real(dp), Intent(IN) :: ucell(3,3) ! Lattice unit cell 696 Character(5), Intent(IN) :: label ! For objects name 697 698! Index for various loops 699 Integer :: i, j 700 701! ---Io/Formats---- 702 200 Format (3F30.25) ! For printing matrix rows (lattice vectors) of ucell 703 210 Format (A, ' ---> ',3F30.25) 704 220 Format (3F15.10) ! A (x,y,z) point 705 706! LATTICE_VECTORS 707 Write(WRITE_UNIT,*) '# ' // trim(label) // ' LATTICE VECTOR INFO' 708 Write(WRITE_UNIT,*) 'object "lattice_vectors" ' // trim(label) // & 709 '" class array type float rank 1 shape 3 items 3 data follows' 710 711! Writing data (cell vectors) 712 Do i=1, 3 713 Write(WRITE_UNIT,200) (ucell(i,j), j=1, 3) 714 End do 715 Write(WRITE_UNIT,*) 'object "lattice_origin' // trim(label) // & 716 '" class array type float rank 1 shape 3 items 3 data follows' 717 718! Writing data (origin of cell) 719 Do i=1, 3 720 Write(WRITE_UNIT,220) 0.0,0.0,0.0 721 End do 722 723 Write(WRITE_UNIT,*) 'object "vectors' // trim(label) // '" class field' 724 Write(WRITE_UNIT,*) 'component "data" value "lattice_vectors' // trim(label) // '"' 725 Write(WRITE_UNIT,*) 'component "positions" value "lattice_origin' // trim(label) // '"' 726 727! UNIT CELL FRAME 728 Write(WRITE_UNIT,*) '# ' // trim(label) // ' UNIT CELL FRAME INFO:' 729 730 Write(WRITE_UNIT,*) 'object "cell_bonds' // trim(label) // '" class array type int rank 1 shape 2 items 12 data follows' 731 732! Writing data [One connection between #a and #b] 733! (N.B: See positions of each corner point below...) 734 Write(WRITE_UNIT,*) ' 0 1' 735 Write(WRITE_UNIT,*) ' 0 2' 736 Write(WRITE_UNIT,*) ' 0 3' 737 Write(WRITE_UNIT,*) ' 1 4' 738 Write(WRITE_UNIT,*) ' 1 5' 739 Write(WRITE_UNIT,*) ' 3 5' 740 Write(WRITE_UNIT,*) ' 3 6' 741 Write(WRITE_UNIT,*) ' 2 6' 742 Write(WRITE_UNIT,*) ' 2 4' 743 Write(WRITE_UNIT,*) ' 7 5' 744 Write(WRITE_UNIT,*) ' 7 6' 745 Write(WRITE_UNIT,*) ' 7 4' 746 747 Write(WRITE_UNIT,*) 'attribute "element type" string "lines"' 748 Write(WRITE_UNIT,*) 'object "cell_corners' // trim(label) // '" class array type float rank 1 shape 3 items 8 data follows' 749 750! Writing data 751 Write(WRITE_UNIT,220) 0.0 , 0.0 , 0.0 ! #0 Origin 752 Write(WRITE_UNIT,220) ucell(1,1), ucell(1,2) , ucell(1,3) ! #1 x 753 Write(WRITE_UNIT,220) ucell(2,1), ucell(2,2) , ucell(2,3) ! #2 y 754 Write(WRITE_UNIT,220) ucell(3,1), ucell(3,2) , ucell(3,3) ! #3 z 755 Write(WRITE_UNIT,220) (ucell(1,1)+ucell(2,1)), & 756 (ucell(1,2)+ucell(2,2)), (ucell(1,3)+ucell(2,3)) ! #4 x+y 757 Write(WRITE_UNIT,220) (ucell(1,1)+ucell(3,1)), & 758 (ucell(1,2)+ucell(3,2)), (ucell(1,3)+ucell(3,3)) ! #5 x+z 759 Write(WRITE_UNIT,220) (ucell(2,1)+ucell(3,1)), & 760 (ucell(2,2)+ucell(3,2)), (ucell(2,3)+ucell(3,3)) ! #6 y+z 761 Write(WRITE_UNIT,220) & 762 (ucell(1,1)+ucell(2,1)+ucell(3,1)), & 763 (ucell(1,2)+ucell(2,2)+ucell(3,2)), & 764 (ucell(1,3)+ucell(2,3)+ucell(3,3)) ! #7 x+y+z 765 766 Write(WRITE_UNIT,*) 'object "bonds' // trim(label) // & 767 '" array type float rank 0 items 12 data follows' 768 769! Writing data 770 Do i=1, 12 771 Write(WRITE_UNIT,*) '1.0' 772 End do 773 Write(WRITE_UNIT,*) 'attribute "dep" string "connections"' 774 775 Write(WRITE_UNIT,*) 'object "cell' // trim(label) // '" class field' 776 Write(WRITE_UNIT,*) 'component "data" value "bonds' // trim(label) // '"' 777 Write(WRITE_UNIT,*) 'component "positions" value "cell_corners' // trim(label) // '"' 778 Write(WRITE_UNIT,*) 'component "connections" value "cell_bonds' // trim(label) // '"' 779 Write(WRITE_UNIT,*) 780 781 End subroutine 782 783! ------------------------------------------------------------------ 784! Subroutines for printing usage and program description in terminal 785! ------------------------------------------------------------------ 786 Subroutine printUsage () 787 Implicit none 788 789 Write(*,*) 'USAGE: Dxformat system_label [-s, -t, -d diff_label]' 790 Write(*,*) ' system_label: the one defined in the .fdf input ' 791 Write(*,*) ' used by Siesta for your simulation' 792 Write(*,*) ' Options: -s --> shifting realspace grid points for isosurfaces (1)' 793 Write(*,*) ' -t --> translation of atoms coordinates (2)' 794 Write(*,*) ' -d --> Use files from second system_label' 795 Write(*,*) " 'diff_label' as reference base data (3)" 796 Write(*,*) '(Note: Use DxFormat alone to see program description)' 797 798 End subroutine 799! ------------------------------------------------------------------ 800 Subroutine printInfo () 801 Implicit none 802 803 Write(*,*) 'This program converts Siesta output from .RHO, .DRHO, ' 804 Write(*,*) '.VH, .VT, .IOCH, .TOCH, .LDOS and .XV files in .dx ' 805 Write(*,*) '.dx format and combine them in one big file for use ' 806 Write(*,*) 'with OpenDx visualization program DxView.net ' 807 Write(*,*) 808 Write(*,*) 'Notes on options:' 809 Write(*,*) ' 1) Siesta uses periodic cell representation and some ' 810 Write(*,*) ' shifting corrections must in general be done for' 811 Write(*,*) ' isosurfaces to ensure centered visualization ' 812 Write(*,*) ' 2) Atoms and isosurface are not writen in the same ' 813 Write(*,*) ' coordinates system by Siesta and translation must' 814 Write(*,*) ' be done to ensure consistency with isosurface' 815 Write(*,*) ' (Unfortunately this works only for orthogonal cell' 816 Write(*,*) ' for now...but this will eventually be fixed) ' 817 Write(*,*) ' 3) When comparing your results with your initial ' 818 Write(*,*) ' configuration, viewing only change data can be ' 819 Write(*,*) ' useful, so you can substract from your data any ' 820 Write(*,*) ' configuration that you want to use as a base ' 821 Write(*,*) ' reference to view only the changes' 822 Write(*,*)' (Note: The difference is applied before conversion' 823 Write(*,*) ' to .dx file format for OpenDx)' 824 Write(*,*) 825 Write(*,*) 'Special note: If the file systemLabel.xyz is found, the' 826 Write(*,*) ' legend labels will be replaced by their ' 827 Write(*,*) ' real labels instead of their default id code' 828 Write(*,*) 829 Write(*,*) 'For further info please read the README file' 830 Write(*,*) 831 832 End subroutine 833! ------------------------------------------------------------------ 834 835 836 837