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