Scripts/loadsave.rb

   1 #
   2 #  loadsave.rb
   3 #
   4 #  Created by Toshi Nagata.
   5 #  Copyright 2008 Toshi Nagata. All rights reserved.
   6 #
   7 # This program is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation version 2 of the License.
  10 #
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 require "scanf.rb"
  17
  18 class Molecule
  19
  20   def loadcrd(filename)
  21     if natoms == 0
  22       raise MolbyError, "cannot load crd; the molecule is empty"
  23     end
  24         fp = open(filename, "rb")
  25     count = 0
  26         frame = 0
  27         coords = (0...natoms).collect { Vector3D[0, 0, 0] }
  28     periodic = (self.box && self.box[4].all? { |n| n != 0 })
  29         show_progress_panel("Loading AMBER crd file...")
  30 #    puts "sframe = #{sframe}, pos = #{fp.pos}"
  31     line = fp.gets   #  Skip first line
  32     while 1
  33       line = fp.gets
  34       if line == nil
  35             if (count > 0)
  36                   raise MolbyError, sprintf("crd format error - file ended in the middle of frame %d", frame)
  37                 end
  38             fp.close
  39                 break
  40       end
  41       line.chomp!
  42 #      values = line.split(' ')
  43           values = line.scan(/......../)
  44       if count + values.size > natoms * 3
  45         raise MolbyError, sprintf("crd format error - too many values at line %d in file %s; number of atoms = %d, current frame = %d", fp.lineno, fp.path, natoms, frame)
  46       end
  47           values.each { |v|
  48             coords[count / 3][count % 3] = Float(v)
  49             count += 1
  50           }
  51       if count == natoms * 3
  52             #  End of frame
  53                 if frame == 0 && atoms.all? { |ap| ap.r.x == 0.0 && ap.r.y == 0.0 && ap.r.z == 0.0 }
  54                         #  Do not create a new frame
  55                         atoms.each_with_index { |ap, i| ap.r = coords[i] }
  56                 else
  57                         create_frame([coords])
  58                 end
  59                 if periodic
  60                   #  Should have box information
  61                   line = fp.gets
  62                   if line == nil || (values = line.chomp.split(' ')).length != 3
  63                     raise "The molecule has a periodic cell but the crd file does not contain cell information"
  64               end
  65                   self.cell = [Float(values[0]), Float(values[1]), Float(values[2]), 90, 90, 90]
  66                 end
  67         count = 0
  68         frame += 1
  69                 if frame % 5 == 0
  70                   set_progress_message("Loading AMBER crd file...\n(#{frame} frames completed)")
  71                 end
  72       end
  73     end
  74         hide_progress_panel
  75         if frame > 0
  76           self.frame = self.nframes - 1
  77           return true
  78         else
  79           return false
  80         end
  81   end
  82
  83   def savecrd(filename)
  84     if natoms == 0
  85       raise MolbyError, "cannot save crd; the molecule is empty"
  86     end
  87     fp = open(filename, "wb")
  88         show_progress_panel("Saving AMBER crd file...")
  89         fp.printf("TITLE: %d atoms\n", natoms)
  90         cframe = self.frame
  91         nframes = self.nframes
  92         j = 0
  93         self.update_enabled = false
  94         begin
  95                 (0...nframes).each { |i|
  96                   select_frame(i)
  97                   j = 0
  98                   while (j < natoms * 3)
  99                         w = atoms[j / 3].r[j % 3]
 100                         fp.printf(" %7.3f", w)
 101                         fp.print("\n") if (j % 10 == 9 || j == natoms * 3 - 1)
 102                         j += 1
 103                   end
 104                   if i % 5 == 0
 105                         set_progress_message("Saving AMBER crd file...\n(#{i} frames completed)")
 106                   end
 107                 }
 108         ensure
 109                 self.update_enabled = true
 110         end
 111         select_frame(cframe)
 112         fp.close
 113         hide_progress_panel
 114         true
 115   end
 116
 117   alias :loadmdcrd :loadcrd
 118   alias :savemdcrd :savecrd
 119
 120   def loadlog(filename)
 121
 122     if natoms == 0
 123                 new_unit = true
 124         else
 125                 new_unit = false
 126     end
 127 #       save_undo_enabled = self.undo_enabled?
 128 #       self.undo_enabled = false
 129         self.update_enabled = false
 130         mes = "Loading GAMESS log file"
 131         show_progress_panel(mes)
 132         ne_alpha = ne_beta = 0   #  number of electrons
 133         rflag = nil  #  0, UHF; 1, RHF; 2, ROHF
 134         mo_count = 0
 135         ncomps = 0   #  Number of AO terms per one MO (sum of the number of components over all shells)
 136         alpha_beta = nil   #  Flag to read alpha/beta MO appropriately
 137         begin
 138                 fp = open(filename, "rb")
 139                 if nframes > 0
 140                         create_frame
 141                         frame = nframes - 1
 142                 end
 143                 n = 0
 144                 while 1
 145                         line = fp.gets
 146                         if line == nil
 147                                 fp.close
 148                                 break
 149                         end
 150                         line.chomp!
 151                         if line =~ /ATOM\s+ATOMIC\s+COORDINATES/ || line =~ /COORDINATES OF ALL ATOMS ARE/
 152                                 set_progress_message(mes + "\nReading atomic coordinates...")
 153                                 first_line = (line =~ /ATOMIC/)
 154                                 line = fp.gets    #  Skip one line
 155                                 n = 0
 156                                 coords = []
 157                                 names = []
 158                                 while (line = fp.gets) != nil
 159                                         break if line =~ /^\s*$/ || line =~ /END OF ONE/
 160                                         next if line =~ /-----/
 161                                         name, charge, x, y, z = line.split
 162                                         v = Vector3D[x, y, z]
 163                                         coords.push(v * (first_line ? 0.529177 : 1.0))  #  Bohr to angstrom
 164                                         names.push([name, charge])
 165                                 end
 166                                 if new_unit
 167                                         #  Build a new molecule
 168                                         names.each_index { |i|
 169                                                 ap = add_atom(names[i][0])
 170                                                 ap.atomic_number = names[i][1].to_i
 171                                                 ap.atom_type = ap.element
 172                                                 ap.r = coords[i]
 173                                         }
 174                                         #  Find bonds
 175                                         guess_bonds
 176                                 #       atoms.each { |ap|
 177                                 #               j = ap.index
 178                                 #               (j + 1 ... natoms).each { |k|
 179                                 #                       if calc_bond(j, k) < 1.7
 180                                 #                               create_bond(j, k)
 181                                 #                       end
 182                                 #               }
 183                                 #       }
 184                                         new_unit = false
 185                                         create_frame
 186                                 else
 187                                         create_frame([coords])  #  Should not be (coords)
 188                                 end
 189                         elsif line =~ /EQUILIBRIUM GEOMETRY LOCATED/i
 190                                 set_progress_message(mes + "\nReading optimized coordinates...")
 191                                 fp.gets; fp.gets; fp.gets
 192                                 n = 0
 193                                 while (line = fp.gets) != nil
 194                                         break if line =~ /^\s*$/
 195                                         line.chomp
 196                                         atom, an, x, y, z = line.split
 197                                         ap = atoms[n]
 198                                         ap.r = Vector3D[x, y, z]
 199                                         n += 1
 200                                         break if n >= natoms
 201                                 end
 202                                 if ne_alpha > 0 && ne_beta > 0
 203                                         #  Allocate basis set record again, to update the atomic coordinates
 204                                         allocate_basis_set_record(rflag, ne_alpha, ne_beta)
 205                                 end
 206                         elsif line =~ /ATOMIC BASIS SET/
 207                                 while (line = fp.gets)
 208                                         break if line =~ /SHELL\s+TYPE\s+PRIMITIVE/
 209                                 end
 210                                 line = fp.gets
 211                                 i = -1
 212                                 nprims = 0
 213                                 sym = -10  #  undefined
 214                                 ncomps = 0
 215                                 while (line = fp.gets)
 216                                         break if line =~ /TOTAL NUMBER OF BASIS SET/
 217                                         line.chomp!
 218                                         if line =~ /^\s*$/
 219                                           #  End of one shell
 220                                           add_gaussian_orbital_shell(sym, nprims, i)
 221                                           # puts "add_gaussian_orbital_shell #{sym}, #{nprims}, #{i}"
 222                                           nprims = 0
 223                                           sym = -10
 224                                           next
 225                                         end
 226                                         a = line.split
 227                                         if a.length == 1
 228                                           i += 1
 229                                           line = fp.gets  #  Skip the blank line
 230                                           next
 231                                         elsif a.length == 5 || a.length == 6
 232                                           if sym == -10
 233                                                 case a[1]
 234                                                 when "S"
 235                                                   sym = 0; n = 1
 236                                                 when "P"
 237                                                   sym = 1; n = 3
 238                                                 when "L"
 239                                                   sym = -1; n = 4
 240                                                 when "D"
 241                                                   sym = 2; n = 6
 242                                                 else
 243                                                   raise MolbyError, "Unknown gaussian shell type at line #{fp.lineno}"
 244                                                 end
 245                                                 ncomps += n
 246                                           end
 247                                           if (a.length == 5 && sym == -1) || (a.length == 6 && sym != -1)
 248                                             raise MolbyError, "Wrong format in gaussian shell information at line #{fp.lineno}"
 249                                           end
 250                                           exp = Float(a[3])
 251                                           c = Float(a[4])
 252                                           csp = Float(a[5] || 0.0)
 253                                           add_gaussian_primitive_coefficients(exp, c, csp)
 254                                           nprims += 1
 255                                           # puts "add_gaussian_primitive_coefficients #{exp}, #{c}, #{csp}"
 256                                         else
 257                                           raise MolbyError, "Error in reading basis set information at line #{fp.lineno}"
 258                                         end
 259                                 end
 260                         elsif line =~ /NUMBER OF OCCUPIED ORBITALS/
 261                                 line =~ /=\s*(\d+)/
 262                                 n = Integer($1)
 263                                 if line =~ /ALPHA/
 264                                         ne_alpha = n
 265                                 else
 266                                         ne_beta = n
 267                                 end
 268                         elsif line =~ /SCFTYP=(\w+)/
 269                                 scftyp = $1
 270                                 if ne_alpha > 0 && ne_beta > 0
 271                                         rflag = 0
 272                                         case scftyp
 273                                         when "RHF"
 274                                                 rflag = 1
 275                                         when "ROHF"
 276                                                 rflag = 2
 277                                         end
 278                                 end
 279                         elsif line =~ /(ALPHA|BETA)\s*SET/
 280                                 alpha_beta = $1
 281                         elsif line =~ /^\s*(EIGENVECTORS|MOLECULAR ORBITALS)\s*$/
 282                                 if mo_count == 0
 283                                         allocate_basis_set_record(rflag, ne_alpha, ne_beta)
 284                                         # puts "allocate_basis_set_record  #{rflag}, #{ne_alpha}, #{ne_beta}"
 285                                 end
 286                                 mo_count += 1
 287                                 idx = 0
 288                                 line = fp.gets; line = fp.gets;
 289                                 set_progress_message(mes + "\nReading MO coefficients...")
 290                                 while (line = fp.gets) != nil
 291                                         break unless line =~ /^\s*\d/
 292                                         mo_labels = line.split       #  MO numbers (1-based)
 293                                         mo_energies = fp.gets.split
 294                                         mo_symmetries = fp.gets.split
 295                                         mo = mo_labels.map { [] }    #  array of *independent* empty arrays
 296                                         while (line = fp.gets) != nil
 297                                                 break unless line =~ /^\s*\d/
 298                                                 line[14..-1].split.each_with_index { |s, i|
 299                                                         mo[i].push(Float(s))
 300                                                 }
 301                                         end
 302                                         mo.each_with_index { |m, i|
 303                                                 idx = Integer(mo_labels[i]) - 1
 304                                                 set_mo_coefficients(idx + (alpha_beta == "BETA" ? ncomps : 0), Float(mo_energies[i]), m)
 305                                         #       if mo_labels[i] % 8 == 1
 306                                         #               puts "set_mo_coefficients #{idx}, #{mo_energies[i]}, [#{m[0]}, ..., #{m[-1]}]"
 307                                         #       end
 308                                         }
 309                                         if line =~ /^\s*$/ && idx < ncomps - 1
 310                                                 next
 311                                         else
 312                                                 break
 313                                         end
 314                                 end
 315                                 set_progress_message(mes)
 316                         end
 317                 end
 318 #       ensure
 319 #               self.undo_enabled = save_undo_enabled
 320 #        hide_progress_panel
 321 #               self.update_enabled = true
 322         end
 323         if nframes > 0
 324           select_frame(nframes - 1)
 325         end
 326         hide_progress_panel
 327         self.update_enabled = true
 328         (n > 0 ? true : false)
 329   end
 330
 331   def loadxyz(filename)
 332 #       save_undo_enabled = self.undo_enabled?
 333 #       self.undo_enabled = false
 334         fp = open(filename, "rb")
 335         n = 0
 336         coords = []
 337         names = []
 338         cell = nil
 339         while 1
 340           line = fp.gets
 341           if line == nil
 342                 fp.close
 343                 break
 344           end
 345           line.chomp
 346           toks = line.split
 347           if coords.length == 0
 348             #  Allow "number of atoms" line or "number of atoms and crystallographic cell parameter"
 349                 #  (Chem3D xyz format)
 350                 next if toks.length == 1
 351                 if toks.length == 7
 352                   cell = toks[1..6].map { |s| Float(s.sub(/\(\d+\)/, "")) }  #  Remove (xx) and convert to float
 353                   next
 354                 end
 355           end
 356           name, x, y, z = line.split
 357           next if z == nil
 358           x = Float(x.sub(/\(\d+\)/, ""))
 359           y = Float(y.sub(/\(\d+\)/, ""))
 360           z = Float(z.sub(/\(\d+\)/, ""))
 361           r = Vector3D[x, y, z]
 362           coords.push(r)
 363           names.push(name)
 364           n += 1
 365         end
 366         celltr = nil
 367         if cell
 368           self.cell = cell
 369           celltr = self.cell_transform
 370         end
 371         names.each_index { |i|
 372           ap = add_atom(names[i])
 373           names[i] =~ /^([A-Za-z]{1,2})/
 374           element = $1.capitalize
 375           ap.element = element
 376           ap.atom_type = element
 377           if celltr
 378             ap.r = celltr * coords[i]
 379           else
 380             ap.r = coords[i]
 381           end
 382         }
 383         guess_bonds
 384         #  Find bonds
 385 #       atoms.each { |ap|
 386 #         j = ap.index
 387 #         (j + 1 ... natoms).each { |k|
 388 #               if calc_bond(j, k) < 1.7
 389 #               #  create_bond(j, k)
 390 #               end
 391 #         }
 392 #       }
 393 #       self.undo_enabled = save_undo_enabled
 394         (n > 0 ? true : false)
 395   end
 396
 397   def loadout(filename)
 398
 399     if natoms == 0
 400                 new_unit = true
 401  #     raise MolbyError, "cannot load crd; the molecule is empty"
 402         else
 403                 new_unit = false
 404     end
 405 #       save_undo_enabled = undo_enabled?
 406 #       self.undo_enabled = false
 407         fp = open(filename, "rb")
 408         if nframes > 0
 409                 create_frame
 410                 frame = nframes - 1
 411         end
 412         n = 0
 413         nf = 0
 414         use_input_orientation = false
 415         show_progress_panel("Loading Gaussian out file...")
 416         while 1
 417                 line = fp.gets
 418                 if line == nil
 419                         fp.close
 420                         break
 421                 end
 422                 line.chomp
 423                 if line =~ /(Input|Standard) orientation/
 424                         match = $1
 425                         if match == "Input"
 426                                 use_input_orientation = true if nf == 0
 427                                 next if !use_input_orientation
 428                         else
 429                                 next if use_input_orientation
 430                         end
 431                         4.times { line = fp.gets }    #  Skip four lines
 432                         n = 0
 433                         coords = []
 434                         anums = []
 435                         while (line = fp.gets) != nil
 436                                 break if line =~ /-----/
 437                                 num, charge, type, x, y, z = line.split
 438                                 coords.push(Vector3D[x, y, z])
 439                                 anums.push(charge)
 440                                 n += 1
 441                         end
 442                         if new_unit
 443                                 #  Build a new molecule
 444                                 anums.each_index { |i|
 445                                         ap = add_atom("X")
 446                                         ap.atomic_number = anums[i]
 447                                         ap.atom_type = ap.element
 448                                         ap.name = sprintf("%s%d", ap.element, i)
 449                                         ap.r = coords[i]
 450                                 }
 451                                 #  Find bonds
 452                         #       atoms.each { |ap|
 453                         #               j = ap.index
 454                         #               (j + 1 ... natoms).each { |k|
 455                         #                       if calc_bond(j, k) < 1.7
 456                         #                               create_bond(j, k)
 457                         #                       end
 458                         #               }
 459                         #       }
 460                                 guess_bonds
 461                                 new_unit = false
 462                                 create_frame
 463                         else
 464                                 create_frame([coords])  #  Should not be (coords)
 465                         end
 466                         nf += 1
 467                 end
 468         end
 469         hide_progress_panel
 470 #       self.undo_enabled = save_undo_enabled
 471         (n > 0 ? true : false)
 472   end
 473
 474   def loadcom(filename)
 475 #       save_undo_enabled = self.undo_enabled?
 476 #       self.undo_enabled = false
 477         self.remove(All)
 478         fp = open(filename, "rb")
 479         section = 0
 480         while (line = fp.gets)
 481           line.chomp!
 482           if section == 0
 483             section = 1 if line =~ /^\#/
 484                 next
 485           elsif section == 1 || section == 2
 486             section += 1 if line =~ /^\s*$/
 487                 next
 488           else
 489             #  The first line is skipped (charge and multiplicity)
 490                 while (line = fp.gets)
 491                   line.chomp!
 492                   break if line =~ /^\s*$/
 493                   a = line.split(/\s*[ \t,\/]\s*/)
 494                   r = Vector3D[Float(a[1]), Float(a[2]), Float(a[3])]
 495                   ap = add_atom(a[0])
 496                   a[0] =~ /^([A-Za-z]{1,2})/
 497                   element = $1.capitalize
 498               ap.element = element
 499               ap.atom_type = element
 500               ap.r = r
 501                 end
 502                 break
 503           end
 504         end
 505         fp.close
 506 #       self.undo_enabled = save_undo_enabled
 507         return true
 508   end
 509
 510   def loadinp(filename)
 511 #       save_undo_enabled = self.undo_enabled?
 512 #       self.undo_enabled = false
 513         self.remove(All)
 514         fp = open(filename, "rb")
 515         section = 0
 516         has_basis = false
 517         while (line = fp.gets)
 518           if line =~ /\A \$BASIS/
 519             has_basis = true
 520                 next
 521           end
 522           next if line !~ /\A \$DATA/
 523           line = fp.gets   #  Title line
 524           line = fp.gets   #  Symmetry line
 525           while (line = fp.gets)
 526 #           puts line
 527             line.chomp!
 528                 break if line =~ /\$END/
 529                 a = line.split
 530                 ap = add_atom(a[0])
 531                 ap.atomic_number = Integer(a[1])
 532                 ap.atom_type = ap.element
 533                 r = Vector3D[Float(a[2]), Float(a[3]), Float(a[4])]
 534                 ap.r = r;
 535                 if !has_basis
 536                   #  Skip until one blank line is detected
 537                   while (line = fp.gets) && line =~ /\S/
 538                   end
 539                 end
 540       end
 541           break
 542         end
 543         fp.close
 544         guess_bonds
 545 #       self.undo_enabled = save_undo_enabled
 546         return true
 547   end
 548
 549   def saveinp(filename)
 550     if natoms == 0
 551       raise MolbyError, "cannot save GAMESS input; the molecule is empty"
 552     end
 553     fp = open(filename, "wb")
 554         now = Time.now.to_s
 555         fp.print <<end_of_header
 556 !  GAMESS input
 557 !  Generated by Molby at #{now}
 558  $CONTRL COORD=UNIQUE EXETYP=RUN ICHARG=0
 559          ICUT=20 INTTYP=HONDO ITOL=30
 560          MAXIT=200 MOLPLT=.T. MPLEVL=0
 561          MULT=1 QMTTOL=1e-08 RUNTYP=OPTIMIZE
 562          SCFTYP=RHF UNITS=ANGS                  $END
 563  $SCF    CONV=1.0E-06 DIRSCF=.T.                $END
 564  $STATPT NSTEP=400 OPTTOL=1.0E-06               $END
 565  $SYSTEM MEMDDI=0 MWORDS=16 TIMLIM=50000        $END
 566  $BASIS  GBASIS=N31 NDFUNC=1 NGAUSS=6           $END
 567  $GUESS  GUESS=HUCKEL                           $END
 568 !
 569  $DATA
 570  #{name}
 571  C1
 572 end_of_header
 573         each_atom { |ap|
 574                 fp.printf " %-6s %4d %10.6f %10.6f %10.6f\n", ap.name, ap.atomic_number, ap.r.x, ap.r.y, ap.r.z
 575         }
 576         fp.print " $END\n"
 577         fp.close
 578         return true
 579   end
 580
 581   def savecom(filename)
 582     if natoms == 0
 583       raise MolbyError, "cannot save Gaussian input; the molecule is empty"
 584     end
 585     fp = open(filename, "wb")
 586         base = File.basename(filename, ".*")
 587         fp.print <<end_of_header
 588 %Chk=#{base}.chk
 589 \# PM3 Opt
 590
 591  #{name}; created by Molby at #{Time.now.to_s}
 592
 593  0 1
 594 end_of_header
 595         each_atom { |ap|
 596                 fp.printf "%-6s %10.6f %10.6f %10.6f\n", ap.element, ap.r.x, ap.r.y, ap.r.z
 597         }
 598         fp.print "\n"
 599         fp.close
 600         return true
 601   end
 602
 603   alias :loadgjf :loadcom
 604   alias :savegjf :savecom
 605
 606   def loadcif(filename)
 607     def getciftoken(fp)
 608           while @tokens.length == 0
 609             line = fp.gets
 610             return nil if !line
 611                 if line[0] == ?;
 612                   s = line
 613                   while line = fp.gets
 614                     break if line[0] == ?;
 615                     s += line
 616                   end
 617                   return s if !line
 618                   s += ";"
 619                   @tokens.push(s)
 620                   line = line[1...line.length]
 621                 end
 622             line.strip!
 623             line.scan(/'[^\']*'|"[^\"]*"|[^#\s]+|#.*/) do |s|
 624               next if s[0] == ?#
 625                   if s =~ /^data_|loop_|global_|save_|stop_/i
 626                     s = "#" + s    #  Label for reserved words
 627                   end
 628                   @tokens.push(s)
 629             end
 630           end
 631           return @tokens.shift
 632         end
 633         def float_strip_rms(str)
 634           return Float(str.sub(/\(\d+\)/, ""))
 635         end
 636         @tokens = []
 637         self.remove(All)
 638         fp = open(filename, "rb")
 639         cell = []
 640         token = getciftoken(fp)
 641         pardigits_re = /\(\d+\)/
 642         while token != nil
 643           if token =~ /^_cell/
 644                 val = getciftoken(fp)
 645                 if token == "_cell_length_a"
 646                   cell[0] = float_strip_rms(val)
 647                 elsif token == "_cell_length_b"
 648                   cell[1] = float_strip_rms(val)
 649                 elsif token == "_cell_length_c"
 650                   cell[2] = float_strip_rms(val)
 651                 elsif token == "_cell_angle_alpha"
 652                   cell[3] = float_strip_rms(val)
 653                 elsif token == "_cell_angle_beta"
 654                   cell[4] = float_strip_rms(val)
 655                 elsif token == "_cell_angle_gamma"
 656                   cell[5] = float_strip_rms(val)
 657                 end
 658                 if cell.length == 6 && cell.all?
 659                   self.cell = cell
 660                   puts "Unit cell is set to #{cell.inspect}."
 661                   cell = []
 662                 end
 663                 token = getciftoken(fp)
 664                 next
 665       elsif token.casecmp("#loop_") == 0
 666             labels = []
 667                 while (token = getciftoken(fp)) && token[0] == ?_
 668                   labels.push(token)
 669                 end
 670                 if labels[0] =~ /symmetry_equiv_pos|atom_site_label|atom_site_aniso_label|geom_bond/
 671                   hlabel = Hash.new(-10000000)
 672                   labels.each_with_index { |lb, i|
 673                         hlabel[lb] = i
 674                   }
 675                   data = []
 676                   n = labels.length
 677                   a = []
 678                   while 1
 679                         break if token == nil || token[0] == ?_ || token[0] == ?#
 680                         a.push(token)
 681                         if a.length == n
 682                           data.push(a)
 683                           a = []
 684                         end
 685                         token = getciftoken(fp)
 686                   end
 687                   if labels[0] =~ /^_symmetry_equiv_pos/
 688                     data.each { |d|
 689                           symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]]
 690                           symstr.delete("\"\'")
 691                           exps = symstr.split(/,/)
 692                           sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
 693                           exps.each_with_index { |s, i|
 694                             terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
 695                                 terms.each { |a|
 696                                   #  a[0]: sign, a[2]: numerator, a[4]: denometer
 697                                   if a[4] != nil
 698                                     #  The number part is a[2]/a[4]
 699                                     num = Float(a[2])/Float(a[4])
 700                                   elsif a[2] != nil
 701                                     #  The number part is either integer or a floating point
 702                                     num = Float(a[2])
 703                                   else
 704                                     num = 1.0
 705                                   end
 706                                   num = -num if a[0][0] == ?-
 707                                   xyz = (a[5] || a[6])
 708                                   if xyz == "x" || xyz == "X"
 709                                     sym[i * 3] = num
 710                                   elsif xyz == "y" || xyz == "Y"
 711                                     sym[i * 3 + 1] = num
 712                                   elsif xyz == "z" || xyz == "Z"
 713                                     sym[i * 3 + 2] = num
 714                                   else
 715                                     sym[9 + i] = num
 716                                   end
 717                                 }
 718                           }
 719                           puts "symmetry operation #{sym.inspect}"
 720                           add_symmetry(Transform.new(sym))
 721                         }
 722                         puts "#{self.nsymmetries} symmetry operations are added"
 723                   elsif labels[0] =~ /^_atom_site_label/
 724                         #  Create atoms
 725                         data.each { |d|
 726                           name = d[hlabel["_atom_site_label"]]
 727                           elem = d[hlabel["_atom_site_type_symbol"]]
 728                           fx = d[hlabel["_atom_site_fract_x"]]
 729                           fy = d[hlabel["_atom_site_fract_y"]]
 730                           fz = d[hlabel["_atom_site_fract_z"]]
 731                           uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
 732                           biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
 733                           occ = d[hlabel["_atom_site_occupancy"]]
 734                           ap = self.add_atom(name, elem, elem)
 735                           ap.fract_x = float_strip_rms(fx)
 736                           ap.fract_y = float_strip_rms(fy)
 737                           ap.fract_z = float_strip_rms(fz)
 738                           if biso
 739                             ap.temp_factor = float_strip_rms(biso)
 740                           elsif uiso
 741                             ap.temp_factor = float_strip_rms(uiso) * 78.9568352087149 #  8*pi*pi
 742                           end
 743                           ap.occupancy = float_strip_rms(occ)
 744                         }
 745                         puts "#{self.natoms} atoms are created."
 746                   elsif labels[0] =~ /^_atom_site_aniso_label/
 747                     #  Set anisotropic parameters
 748                         c = 0
 749                         data.each { |d|
 750                           name = d[hlabel["_atom_site_aniso_label"]]
 751                           ap = self.atoms[name]
 752                           next if !ap
 753                           u11 = d[hlabel["_atom_site_aniso_U_11"]]
 754                           if u11
 755                             u11 = float_strip_rms(u11)
 756                             u22 = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
 757                             u33 = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
 758                             u12 = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
 759                             u13 = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
 760                             u23 = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
 761                             ap.aniso = [u11, u22, u33, u12, u13, u23, 8]
 762                                 c += 1
 763                           end
 764                         }
 765                         puts "#{c} anisotropic parameters are set."
 766                   elsif labels[0] =~ /^_geom_bond/
 767                     #  Create bonds
 768                         exbonds = []
 769                         data.each { |d|
 770                           n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
 771                           n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
 772                           sym1 = d[hlabel["_geom_bond_site_symmetry_1"]]
 773                           sym2 = d[hlabel["_geom_bond_site_symmetry_2"]]
 774                           if sym1 != "." || sym2 != "."
 775                             exbonds.push([n1, n2, sym1, sym2])
 776                           else
 777                             self.create_bond(n1, n2)
 778                           end
 779                     }
 780                         if exbonds.length > 0
 781                           h = Dialog.run {
 782                             layout(1,
 783                                   item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
 784                                   item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
 785                                   item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
 786                                   item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
 787                                 )
 788                                 radio_group("atoms_only", "fragment", "ignore")
 789                           }
 790                           if h[:status] == 0 && h["ignore"] == 0
 791                             atoms_only = (h["atoms_only"] != 0)
 792                                 if !atoms_only
 793                                   fragments = []
 794                                   self.each_fragment { |f| fragments.push(f) }
 795                                 end
 796                                 sym_decoder = /(\d+)_(\d)(\d)(\d)/
 797                                 debug = nil
 798                                 exbonds.each { |ex|
 799                                   #  Convert name to index before expansion
 800                                   if debug
 801                                     ex[4] = ex[0]
 802                                         ex[5] = ex[1]
 803                                   end
 804                                   ex[0] = self.atoms[ex[0]].index
 805                                   ex[1] = self.atoms[ex[1]].index
 806                                 }
 807                                 exbonds.each { |ex|
 808                                   if debug; puts "extra bond #{ex[4]}(#{ex[2]}) - #{ex[5]}(#{ex[3]})"; end
 809                                   (2..3).each { |i|
 810                                     if ex[i] == "."
 811                                           ex[i] = ex[i - 2]    #  No expansion
 812                                         elsif ex[i] =~ /(\d+)_(\d)(\d)(\d)/
 813                                   symop = [Integer($1) - 1, Integer($2) - 5, Integer($3) - 5, Integer($4) - 5, ex[i - 2]]
 814                                           if debug; puts "  symop = #{symop.inspect}"; end
 815                                           ap = self.atoms.find { |ap| (s = ap.symop) != nil && s === symop }
 816                                           if ap
 817                                             if debug; puts "  already expanded (atom #{ap.index}, #{ap.name}, #{ap.symop.inspect})"; end
 818                                             ex[i] = ap.index   #  Already expanded
 819                                           else
 820                                             #  Expand the atom or the fragment including the atom
 821                                                 if atoms_only
 822                                                   ig = IntGroup[ex[i - 2]]
 823                                                 else
 824                                                   ig = fragments.find { |f| f.include?(ex[i - 2]) }
 825                                                 end
 826                                             if debug; puts "  expanding #{ig} by #{symop.inspect}"; end
 827                                                 self.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
 828                                                 #  Find again the expanded atom
 829                                             ap = self.atoms.find { |ap| (s = ap.symop) != nil && s === symop }
 830                                                 ex[i] = ap.index
 831                                           end
 832                                         else
 833                                           raise "unrecognizable symmetry operation: #{ex[i]}"
 834                                         end
 835                                   }
 836                                   if debug; puts "  creating bond #{ex[2]} - #{ex[3]}"; end
 837                                   self.create_bond(ex[2], ex[3])
 838                                 }
 839                           end
 840                         end
 841                         puts "#{self.nbonds} bonds are created."
 842                   end
 843                   next
 844                 else
 845                 #  puts "Loop beginning with #{labels[0]} is skipped"
 846                 end
 847           else
 848             #  Skip this token
 849                 token = getciftoken(fp)
 850           end
 851           #  Skip tokens until next tag or reserved word is detected
 852           while token != nil && token[0] != ?_ && token[0] != ?#
 853                 token = getciftoken(fp)
 854           end
 855           next
 856         end
 857         fp.close
 858 #       self.undo_enabled = save_undo_enabled
 859         return true
 860   end
 861
 862   def dump(group = nil)
 863     def quote(str)
 864           if str == ""
 865             str = "\"\""
 866           else
 867             str = str.gsub(/%/, "%%")
 868                 str.gsub!(/ /, "%20")
 869                 str.gsub!(/\t/, "%09")
 870           end
 871           str
 872         end
 873         group = atom_group(group ? group : 0...natoms)
 874         s = ""
 875         group.each { |i|
 876           ap = atoms[i]
 877           s += sprintf("%4d %-7s %-4s %-4s %-2s %7.3f %7.3f %7.3f %6.3f [%s]\n",
 878                 ap.index, sprintf("%3s.%d", ap.res_name, ap.res_seq),
 879                 quote(ap.name), quote(ap.atom_type), quote(ap.element),
 880                 ap.r.x, ap.r.y, ap.r.z, ap.charge,
 881                 ap.connects.join(","))
 882         }
 883         print s
 884   end
 885
 886   def from_dump(arg)
 887     if natoms > 0
 888           raise "Molecule must be empty"
 889         end
 890     format = "index residue name atom_type element rx ry rz charge connects"
 891         keys = []
 892         resAtoms = Hash.new
 893         newBonds = []
 894     arg.each { |line|
 895       #  arg can be either a String or an array of String. If it is a string,
 896           #  arg.each iterates for each line in the string. If it is an array,
 897           #  arg.each iterates for each member of the array.
 898           if line =~ /^\#/
 899             format = line[1..-1]
 900                 keys = []
 901           end
 902           if keys.length == 0
 903             keys = format.split(" ").collect { |key| key.to_sym }
 904           end
 905           values = line.chomp.split(" ")
 906           next if values == nil || values.length == 0
 907           ap = create_atom(sprintf("X%03d", natoms))
 908           r = Vector3D[0, 0, 0]
 909           keys.each_index { |i|
 910             break if (value = values[i]) == nil
 911                 if value == "\"\""
 912                   value = ""
 913                 else
 914                   value.gsub(/%09/, "\t")
 915                   value.gsub(/%20/, " ")
 916                   value.gsub(/%%/, "%")
 917                 end
 918                 key = keys[i]
 919                 if key == :residue
 920                   if resAtoms[value] == nil
 921                     resAtoms[value] = []
 922                   end
 923                   resAtoms[value].push(ap.index)
 924                 elsif key == :rx
 925                   r.x = value.to_f
 926                 elsif key == :ry
 927                   r.y = value.to_f
 928                 elsif key == :rz
 929                   r.z = value.to_f
 930                 elsif key == :connects
 931                   value.scan(/\d+/).each { |i|
 932                     i = i.to_i
 933                     if ap.index < i
 934                           newBonds.push(ap.index)
 935                           newBonds.push(i)
 936                         end
 937                   }
 938                 elsif key == :index
 939                   next
 940                 else
 941                   ap.set_attr(key, value)
 942                 end
 943           }
 944           ap.r = r
 945         }
 946         resAtoms.each_key { |key|
 947           assign_residue(atom_group(resAtoms[key]), key)
 948         }
 949         create_bond(*newBonds)
 950         self
 951   end
 952
 953 end