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             sym = nil
 678           elsif (str =~ /(\d+)_(\d)(\d)(\d)/) || (str =~ /(\d+) +(\d)(\d)(\d)/)
 679             sym = [Integer($1) - 1, Integer($2) - 5, Integer($3) - 5, Integer($4) - 5]
 680           elsif (str =~ /^(\d+)$/)
 681             sym = [Integer($1) - 1, 0, 0, 0]
 682           end
 683           if sym && (sym[0] == 0 && sym[1] == 0 && sym[2] == 0 && sym[3] == 0)
 684             sym = nil
 685           end
 686           sym
 687         end
 688         warn_message = ""
 689         verbose = nil
 690         bond_defined = false
 691         @tokens = []
 692         special_positions = []
 693         self.remove(All)
 694         fp = open(filename, "rb")
 695         cell = []
 696         cell_trans = cell_trans_inv = Transform.identity
 697         token = getciftoken(fp)
 698         pardigits_re = /\(\d+\)/
 699         calculated_atoms = []
 700         while token != nil
 701           if token =~ /^_cell/
 702                 val = getciftoken(fp)
 703                 if token == "_cell_length_a"
 704                   cell[0], cell[6] = float_strip_rms(val)
 705                 elsif token == "_cell_length_b"
 706                   cell[1], cell[7] = float_strip_rms(val)
 707                 elsif token == "_cell_length_c"
 708                   cell[2], cell[8] = float_strip_rms(val)
 709                 elsif token == "_cell_angle_alpha"
 710                   cell[3], cell[9] = float_strip_rms(val)
 711                 elsif token == "_cell_angle_beta"
 712                   cell[4], cell[10] = float_strip_rms(val)
 713                 elsif token == "_cell_angle_gamma"
 714                   cell[5], cell[11] = float_strip_rms(val)
 715                 end
 716                 if cell.length == 12 && cell.all?
 717                   self.cell = cell
 718                   puts "Unit cell is set to #{cell.inspect}." if verbose
 719                   cell = []
 720                   cell_trans = self.cell_transform
 721                   cell_trans_inv = cell_trans.inverse
 722                 end
 723                 token = getciftoken(fp)
 724                 next
 725       elsif token.casecmp("#loop_") == 0
 726             labels = []
 727                 while (token = getciftoken(fp)) && token[0] == ?_
 728                   labels.push(token)
 729                 end
 730                 if labels[0] =~ /symmetry_equiv_pos|atom_site_label|atom_site_aniso_label|geom_bond/
 731                   hlabel = Hash.new(-10000000)
 732                   labels.each_with_index { |lb, i|
 733                         hlabel[lb] = i
 734                   }
 735                   data = []
 736                   n = labels.length
 737                   a = []
 738                   while 1
 739                         break if token == nil || token[0] == ?_ || token[0] == ?#
 740                         a.push(token)
 741                         if a.length == n
 742                           data.push(a)
 743                           a = []
 744                         end
 745                         token = getciftoken(fp)
 746                   end
 747                   if labels[0] =~ /^_symmetry_equiv_pos/
 748                     data.each { |d|
 749                           symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]]
 750                           symstr.delete("\"\'")
 751                           exps = symstr.split(/,/)
 752                           sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
 753                           exps.each_with_index { |s, i|
 754                             terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
 755                                 terms.each { |a|
 756                                   #  a[0]: sign, a[2]: numerator, a[4]: denometer
 757                                   if a[4] != nil
 758                                     #  The number part is a[2]/a[4]
 759                                     num = Float(a[2])/Float(a[4])
 760                                   elsif a[2] != nil
 761                                     #  The number part is either integer or a floating point
 762                                     num = Float(a[2])
 763                                   else
 764                                     num = 1.0
 765                                   end
 766                                   num = -num if a[0][0] == ?-
 767                                   xyz = (a[5] || a[6])
 768                                   if xyz == "x" || xyz == "X"
 769                                     sym[i] = num
 770                                   elsif xyz == "y" || xyz == "Y"
 771                                     sym[i + 3] = num
 772                                   elsif xyz == "z" || xyz == "Z"
 773                                     sym[i + 6] = num
 774                                   else
 775                                     sym[9 + i] = num
 776                                   end
 777                                 }
 778                           }
 779                           puts "symmetry operation #{sym.inspect}" if verbose
 780                           add_symmetry(Transform.new(sym))
 781                         }
 782                         puts "#{self.nsymmetries} symmetry operations are added" if verbose
 783                   elsif labels[0] =~ /^_atom_site_label/
 784                         #  Create atoms
 785                         data.each { |d|
 786                           name = d[hlabel["_atom_site_label"]]
 787                           elem = d[hlabel["_atom_site_type_symbol"]]
 788                           fx = d[hlabel["_atom_site_fract_x"]]
 789                           fy = d[hlabel["_atom_site_fract_y"]]
 790                           fz = d[hlabel["_atom_site_fract_z"]]
 791                           uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
 792                           biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
 793                           occ = d[hlabel["_atom_site_occupancy"]]
 794                           calc = d[hlabel["_atom_site_calc_flag"]]
 795                           if elem == nil || elem == ""
 796                             if name =~ /[A-Za-z]{1,2}/
 797                                   elem = $&.capitalize
 798                                 else
 799                                   elem = "Du"
 800                                 end
 801                           end
 802                           ap = self.add_atom(name, elem, elem)
 803                           ap.fract_x, ap.sigma_x = float_strip_rms(fx)
 804                           ap.fract_y, ap.sigma_y = float_strip_rms(fy)
 805                           ap.fract_z, ap.sigma_z = float_strip_rms(fz)
 806                           if biso
 807                             ap.temp_factor, sig = float_strip_rms(biso)
 808                           elsif uiso
 809                             ap.temp_factor, sig = float_strip_rms(uiso)
 810                                 ap.temp_factor *= 78.9568352087149          #  8*pi*pi
 811                           end
 812                           ap.occupancy, sig = float_strip_rms(occ)
 813                           if calc == "c" || calc == "calc"
 814                             calculated_atoms.push(ap.index)
 815                       end
 816                           #  Guess special positions
 817                           (1...nsymmetries).each { |isym|
 818                             sr = ap.fract_r
 819                             sr = (transform_for_symop(isym) * sr) - sr;
 820                                 nx = (sr.x + 0.5).floor
 821                                 ny = (sr.y + 0.5).floor
 822                                 nz = (sr.z + 0.5).floor
 823                                 if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
 824                                   #  [isym, -nx, -ny, -nz] transforms this atom to itself
 825                                   #  The following line is equivalent to:
 826                                   #    if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
 827                                   #    special_positions[ap.index].push(...)
 828                                   (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
 829                                 end
 830                           }
 831                           if verbose && special_positions[ap.index]
 832                             puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
 833                           end
 834                         }
 835                         puts "#{self.natoms} atoms are created." if verbose
 836                   elsif labels[0] =~ /^_atom_site_aniso_label/
 837                     #  Set anisotropic parameters
 838                         c = 0
 839                         data.each { |d|
 840                           name = d[hlabel["_atom_site_aniso_label"]]
 841                           ap = self.atoms[name]
 842                           next if !ap
 843                           u11 = d[hlabel["_atom_site_aniso_U_11"]]
 844                           if u11
 845                             usig = []
 846                             u11, usig[0] = float_strip_rms(u11)
 847                             u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
 848                             u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
 849                             u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
 850                             u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
 851                             u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
 852                             ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
 853                                 c += 1
 854                           end
 855                         }
 856                         puts "#{c} anisotropic parameters are set." if verbose
 857                   elsif labels[0] =~ /^_geom_bond/
 858                     #  Create bonds
 859                         exbonds = []
 860                         data.each { |d|
 861                           n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
 862                           n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
 863                           sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
 864                           sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
 865                           n1 = self.atoms[n1].index
 866                           n2 = self.atoms[n2].index
 867                           sym1 = parse_symmetry_operation(sym1)
 868                           sym2 = parse_symmetry_operation(sym2)
 869                           if sym1 || sym2
 870                             exbonds.push([n1, n2, sym1, sym2])
 871                           else
 872                             self.create_bond(n1, n2)
 873                           end
 874                           tr1 = (sym1 ? transform_for_symop(sym1) : Transform.identity)
 875                           tr2 = (sym2 ? transform_for_symop(sym2) : Transform.identity)
 876                           if special_positions[n1]
 877                                 #  Add extra bonds for equivalent positions of n1
 878                                 special_positions[n1].each { |symop|
 879                                   sym2x = symop_for_transform(tr1 * transform_for_symop(symop) * tr1.inverse * tr2)
 880                                   exbonds.push([n1, n2, sym1, sym2x])
 881                                 }
 882                           end
 883                           if special_positions[n2]
 884                                 #  Add extra bonds n2-n1.symop, where symop transforms n2 to self
 885                                 tr = (sym1 ? transform_for_symop(sym1) : Transform.identity)
 886                                 special_positions[n2].each { |symop|
 887                                   sym1x = symop_for_transform(tr2 * transform_for_symop(symop) * tr2.inverse * tr1)
 888                                   exbonds.push([n2, n1, sym2, sym1x])
 889                                 }
 890                           end
 891                     }
 892                         bond_defined = true
 893                         puts "#{self.nbonds} bonds are created." if verbose
 894                         if calculated_atoms.length > 0
 895                           #  Guess bonds for calculated hydrogen atoms
 896                           n1 = 0
 897                           calculated_atoms.each { |ai|
 898                             if atoms[ai].connects.length == 0
 899                                   as = find_close_atoms(ai)
 900                                   as.each { |aj|
 901                                     self.create_bond(ai, aj)
 902                                         n1 += 1
 903                                   }
 904                                 end
 905                           }
 906                           puts "#{n1} bonds are guessed." if verbose
 907                         end
 908                         if exbonds.length > 0
 909                           h = Dialog.run {
 910                             layout(1,
 911                                   item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
 912                                   item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
 913                                   item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
 914                                   item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
 915                                 )
 916                                 radio_group("atoms_only", "fragment", "ignore")
 917                           }
 918                           if h[:status] == 0 && h["ignore"] == 0
 919                             atoms_only = (h["atoms_only"] != 0)
 920                                 if !atoms_only
 921                                   fragments = []
 922                                   self.each_fragment { |f| fragments.push(f) }
 923                                 end
 924                                 debug = nil
 925                                 exbonds.each { |ex|
 926                                   if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
 927                                   ex0 = ex.dup
 928                                   (2..3).each { |i|
 929                                     symop = ex[i]
 930                                         if symop == nil
 931                                           ex[i + 2] = ex[i - 2]
 932                                         else
 933                                           if debug; puts "  symop = #{symop.inspect}"; end
 934                                           #  Expand the atom or the fragment including the atom
 935                                           if atoms_only
 936                                                 ig = IntGroup[ex[i - 2]]
 937                                                 idx = 0
 938                                           else
 939                                                 ig = fragments.find { |f| f.include?(ex[i - 2]) }
 940                                                 ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
 941                                           end
 942                                           symop[4] = ex[i - 2]  #  Base atom
 943                                           if debug; puts "  expanding #{ig} by #{symop.inspect}"; end
 944                                           a = self.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
 945                                           ex[i + 2] = a[idx]   #  Index of the expanded atom
 946                                         end
 947                                   }
 948                                   if ex[4] && ex[5] && ex[4] != ex[5]
 949                                     if debug; puts "  creating bond #{ex[4]} - #{ex[5]}"; end
 950                                     self.create_bond(ex[4], ex[5])
 951                                   end
 952                                 }
 953                           end
 954                         end
 955                         puts "#{self.nbonds} bonds are created." if verbose
 956                   end
 957                   next
 958                 else
 959                 #  puts "Loop beginning with #{labels[0]} is skipped"
 960                 end
 961           else
 962             #  Skip this token
 963                 token = getciftoken(fp)
 964           end
 965           #  Skip tokens until next tag or reserved word is detected
 966           while token != nil && token[0] != ?_ && token[0] != ?#
 967                 token = getciftoken(fp)
 968           end
 969           next
 970         end
 971         fp.close
 972         if !bond_defined
 973           self.guess_bonds
 974         end
 975 #       self.undo_enabled = save_undo_enabled
 976         return true
 977   end
 978
 979   def dump(group = nil)
 980     def quote(str)
 981           if str == ""
 982             str = "\"\""
 983           else
 984             str = str.gsub(/%/, "%%")
 985                 str.gsub!(/ /, "%20")
 986                 str.gsub!(/\t/, "%09")
 987           end
 988           str
 989         end
 990         group = atom_group(group ? group : 0...natoms)
 991         s = ""
 992         group.each { |i|
 993           ap = atoms[i]
 994           s += sprintf("%4d %-7s %-4s %-4s %-2s %7.3f %7.3f %7.3f %6.3f [%s]\n",
 995                 ap.index, sprintf("%3s.%d", ap.res_name, ap.res_seq),
 996                 quote(ap.name), quote(ap.atom_type), quote(ap.element),
 997                 ap.r.x, ap.r.y, ap.r.z, ap.charge,
 998                 ap.connects.join(","))
 999         }
1000         print s
1001   end
1002
1003   def from_dump(arg)
1004     if natoms > 0
1005           raise "Molecule must be empty"
1006         end
1007     format = "index residue name atom_type element rx ry rz charge connects"
1008         keys = []
1009         resAtoms = Hash.new
1010         newBonds = []
1011     arg.each { |line|
1012       #  arg can be either a String or an array of String. If it is a string,
1013           #  arg.each iterates for each line in the string. If it is an array,
1014           #  arg.each iterates for each member of the array.
1015           if line =~ /^\#/
1016             format = line[1..-1]
1017                 keys = []
1018           end
1019           if keys.length == 0
1020             keys = format.split(" ").collect { |key| key.to_sym }
1021           end
1022           values = line.chomp.split(" ")
1023           next if values == nil || values.length == 0
1024           ap = create_atom(sprintf("X%03d", natoms))
1025           r = Vector3D[0, 0, 0]
1026           keys.each_index { |i|
1027             break if (value = values[i]) == nil
1028                 if value == "\"\""
1029                   value = ""
1030                 else
1031                   value.gsub(/%09/, "\t")
1032                   value.gsub(/%20/, " ")
1033                   value.gsub(/%%/, "%")
1034                 end
1035                 key = keys[i]
1036                 if key == :residue
1037                   if resAtoms[value] == nil
1038                     resAtoms[value] = []
1039                   end
1040                   resAtoms[value].push(ap.index)
1041                 elsif key == :rx
1042                   r.x = value.to_f
1043                 elsif key == :ry
1044                   r.y = value.to_f
1045                 elsif key == :rz
1046                   r.z = value.to_f
1047                 elsif key == :connects
1048                   value.scan(/\d+/).each { |i|
1049                     i = i.to_i
1050                     if ap.index < i
1051                           newBonds.push(ap.index)
1052                           newBonds.push(i)
1053                         end
1054                   }
1055                 elsif key == :index
1056                   next
1057                 else
1058                   ap.set_attr(key, value)
1059                 end
1060           }
1061           ap.r = r
1062         }
1063         resAtoms.each_key { |key|
1064           assign_residue(atom_group(resAtoms[key]), key)
1065         }
1066         create_bond(*newBonds)
1067         self
1068   end
1069
1070 end