Scripts/loadsave.rb

   1 # coding: utf-8
   2 #
   3 #  loadsave.rb
   4 #
   5 #  Created by Toshi Nagata.
   6 #  Copyright 2008 Toshi Nagata. All rights reserved.
   7 #
   8 # This program is free software; you can redistribute it and/or modify
   9 # it under the terms of the GNU General Public License as published by
  10 # the Free Software Foundation version 2 of the License.
  11 #
  12 # This program is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16
  17 require "scanf.rb"
  18
  19 class Molecule
  20
  21   def loadcrd(filename)
  22     if natoms == 0
  23       raise MolbyError, "cannot load crd; the molecule is empty"
  24     end
  25         fp = open(filename, "rb")
  26     count = 0
  27         frame = 0
  28         coords = (0...natoms).collect { Vector3D[0, 0, 0] }
  29     periodic = (self.box && self.box[4].any? { |n| n != 0 })
  30         show_progress_panel("Loading AMBER crd file...")
  31 #    puts "sframe = #{sframe}, pos = #{fp.pos}"
  32     while line = fp.gets
  33       line.chomp!
  34           values = line.scan(/......../)
  35           if fp.lineno == 1
  36             #  The first line should be skipped. However, if this line seems to contain
  37                 #  coordinates, then try reading them
  38                 if values.size != 10 && (values.size > 10 || values.size != natoms * 3)
  39                   next  #  Wrong number of coordinates
  40                 end
  41                 if values.each { |v| Float(v) rescue break } == nil
  42                   #  The line contains non-number
  43                   next
  44             end
  45                 puts "Loadcrd: The title line seems to be missing"
  46           end
  47       if count + values.size > natoms * 3
  48         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)
  49       end
  50           values.each { |v|
  51             coords[count / 3][count % 3] = Float(v)
  52             count += 1
  53           }
  54       if count == natoms * 3
  55             #  End of frame
  56                 if frame == 0 && atoms.all? { |ap| ap.r.x == 0.0 && ap.r.y == 0.0 && ap.r.z == 0.0 }
  57                         #  Do not create a new frame
  58                         atoms.each_with_index { |ap, i| ap.r = coords[i] }
  59                 else
  60                         create_frame([coords])
  61                 end
  62                 if periodic
  63                   #  Should have box information
  64                   line = fp.gets
  65                   if line == nil || (values = line.chomp.scan(/......../)).length != 3
  66                     #  Periodic but no cell information
  67                         puts "Loadcrd: the molecule has a periodic cell but the crd file does not contain cell info"
  68                         periodic = false
  69                         redo
  70               end
  71                   self.cell = [Float(values[0]), Float(values[1]), Float(values[2]), 90, 90, 90]
  72                 end
  73         count = 0
  74         frame += 1
  75                 if frame % 5 == 0
  76                   set_progress_message("Loading AMBER crd file...\n(#{frame} frames completed)")
  77                 end
  78       end
  79     end
  80         fp.close
  81         if (count > 0)
  82           raise MolbyError, sprintf("crd format error - file ended in the middle of frame %d", frame)
  83         end
  84         hide_progress_panel
  85         if frame > 0
  86           self.frame = self.nframes - 1
  87           return true
  88         else
  89           return false
  90         end
  91   end
  92
  93   def savecrd(filename)
  94     if natoms == 0
  95       raise MolbyError, "cannot save crd; the molecule is empty"
  96     end
  97     fp = open(filename, "wb")
  98         show_progress_panel("Saving AMBER crd file...")
  99         fp.printf("TITLE: %d atoms\n", natoms)
 100         cframe = self.frame
 101         nframes = self.nframes
 102         j = 0
 103         self.update_enabled = false
 104         begin
 105                 (0...nframes).each { |i|
 106                   select_frame(i)
 107                   j = 0
 108                   while (j < natoms * 3)
 109                         w = atoms[j / 3].r[j % 3]
 110                         fp.printf(" %7.3f", w)
 111                         fp.print("\n") if (j % 10 == 9 || j == natoms * 3 - 1)
 112                         j += 1
 113                   end
 114                   if i % 5 == 0
 115                         set_progress_message("Saving AMBER crd file...\n(#{i} frames completed)")
 116                   end
 117                 }
 118         ensure
 119                 self.update_enabled = true
 120         end
 121         select_frame(cframe)
 122         fp.close
 123         hide_progress_panel
 124         true
 125   end
 126
 127   alias :loadmdcrd :loadcrd
 128   alias :savemdcrd :savecrd
 129
 130   def sub_load_gamess_log_basis_set(lines, lineno)
 131     ln = 0
 132         while (line = lines[ln])
 133                 ln += 1
 134                 break if line =~ /SHELL\s+TYPE\s+PRIMITIVE/
 135         end
 136         ln += 1
 137         i = -1
 138         nprims = 0
 139         sym = -10  #  undefined
 140         ncomps = 0
 141         clear_basis_set
 142         while (line = lines[ln])
 143                 ln += 1
 144                 break if line =~ /TOTAL NUMBER OF BASIS SET/
 145                 if line =~ /^\s*$/
 146                   #  End of one shell
 147                   add_gaussian_orbital_shell(i, sym, nprims)
 148                   nprims = 0
 149                   sym = -10
 150                   next
 151                 end
 152                 a = line.split
 153                 if a.length == 1
 154                   i += 1
 155                   ln += 1  #  Skip the blank line
 156                   next
 157                 elsif a.length == 5 || a.length == 6
 158                   if sym == -10
 159                         case a[1]
 160                         when "S"
 161                           sym = 0; n = 1
 162                         when "P"
 163                           sym = 1; n = 3
 164                         when "L"
 165                           sym = -1; n = 4
 166                         when "D"
 167                           sym = 2; n = 6
 168                         when "F"
 169                           sym = 3; n = 10
 170                         when "G"
 171                           sym = 4; n = 15
 172                         else
 173                           raise MolbyError, "Unknown gaussian shell type at line #{lineno + ln}"
 174                         end
 175                         ncomps += n
 176                   end
 177                   if (a.length == 5 && sym == -1) || (a.length == 6 && sym != -1)
 178                         raise MolbyError, "Wrong format in gaussian shell information at line #{lineno + ln}"
 179                   end
 180                   exp = Float(a[3])
 181                   c = Float(a[4])
 182                   csp = Float(a[5] || 0.0)
 183                   add_gaussian_primitive_coefficients(exp, c, csp)
 184                   nprims += 1
 185                 else
 186                   raise MolbyError, "Error in reading basis set information at line #{lineno + ln}"
 187                 end
 188         end
 189         return ncomps
 190   end
 191
 192   def sub_load_gamess_log_mo_coefficients(lines, lineno, ncomps)
 193     ln = 0
 194         idx = 0
 195         alpha = true
 196         while (line = lines[ln]) != nil
 197                 ln += 1
 198                 if line =~ /BETA SET/
 199                         alpha = false
 200                         next
 201                 end
 202                 if line =~ /------------/ || line =~ /EIGENVECTORS/ || line =~ /\*\*\*\* (ALPHA|BETA) SET/
 203                         next
 204                 end
 205                 next unless line =~ /^\s*\d/
 206                 mo_labels = line.split       #  MO numbers (1-based)
 207                 mo_energies = lines[ln].split
 208                 mo_symmetries = lines[ln + 1].split
 209         #       puts "mo #{mo_labels.inspect}"
 210                 ln += 2
 211                 mo = mo_labels.map { [] }    #  array of *independent* empty arrays
 212                 while (line = lines[ln]) != nil
 213                         ln += 1
 214                         break unless line =~ /^\s*\d/
 215                         5.times { |i|
 216                           s = line[15 + 11 * i, 11].chomp
 217                           break if s =~ /^\s*$/
 218                           mo[i].push(Float(s)) rescue print "line = #{line}, s = #{s}"
 219                         # line[15..-1].split.each_with_index { |s, i|
 220                         #       mo[i].push(Float(s))
 221                         }
 222                 end
 223                 mo.each_with_index { |m, i|
 224                         idx = Integer(mo_labels[i])
 225                         set_mo_coefficients(idx + (alpha ? 0 : ncomps), Float(mo_energies[i]), m)
 226                 #       if mo_labels[i] % 8 == 1
 227                 #               puts "set_mo_coefficients #{idx}, #{mo_energies[i]}, [#{m[0]}, ..., #{m[-1]}]"
 228                 #       end
 229                 }
 230 #               if line =~ /^\s*$/
 231 #                       next
 232 #               else
 233 #                       break
 234 #               end
 235         end
 236   end
 237
 238   def sub_load_gamess_log(fp)
 239
 240     if natoms == 0
 241                 new_unit = true
 242         else
 243                 new_unit = false
 244     end
 245         self.update_enabled = false
 246         mes = "Loading GAMESS log file"
 247         show_progress_panel(mes)
 248         energy = nil
 249         ne_alpha = ne_beta = 0   #  number of electrons
 250         rflag = nil  #  0, UHF; 1, RHF; 2, ROHF
 251         mo_count = 0
 252         search_mode = 0   #  0, no search; 1, optimize; 2, irc
 253         ncomps = 0   #  Number of AO terms per one MO (sum of the number of components over all shells)
 254         alpha_beta = nil   #  Flag to read alpha/beta MO appropriately
 255         nsearch = 0  #  Search number for optimization
 256         begin
 257         #       if nframes > 0
 258         #               create_frame
 259         #               frame = nframes - 1
 260         #       end
 261                 n = 0
 262                 while 1
 263                         line = fp.gets
 264                         if line == nil
 265                                 break
 266                         end
 267                         line.chomp!
 268                         if line =~ /ATOM\s+ATOMIC\s+COORDINATES/ || line =~ /COORDINATES OF ALL ATOMS ARE/ || line =~ /COORDINATES \(IN ANGSTROM\) FOR \$DATA GROUP ARE/
 269                                 set_progress_message(mes + "\nReading atomic coordinates...")
 270                                 if line =~ /ATOMIC/
 271                                   first_line = true
 272                                 #  if !new_unit
 273                                 #    next   #  Skip initial atomic coordinates unless loading into an empty molecule
 274                                 #  end
 275                                   line = fp.gets  #  Skip one line
 276                                 else
 277                                   first_line = false
 278                                   nsearch += 1
 279                                   if line =~ /COORDINATES OF ALL ATOMS ARE/
 280                                     line = fp.gets  #  Skip one line
 281                               end
 282                                 end
 283                                 n = 0
 284                                 coords = []
 285                                 names = []
 286                                 while (line = fp.gets) != nil
 287                                         break if line =~ /^\s*$/ || line =~ /END OF ONE/ || line =~ /\$END/
 288                                         next if line =~ /-----/
 289                                         name, charge, x, y, z = line.split
 290                                         v = Vector3D[x, y, z]
 291                                         coords.push(v * (first_line ? 0.529177 : 1.0))  #  Bohr to angstrom
 292                                         names.push([name, charge])
 293                                 end
 294                                 if new_unit
 295                                         #  Build a new molecule
 296                                         names.each_index { |i|
 297                                                 ap = add_atom(names[i][0])
 298                                                 ap.atomic_number = names[i][1].to_i
 299                                                 ap.atom_type = ap.element
 300                                                 ap.r = coords[i]
 301                                         }
 302                                         #  Find bonds
 303                                         guess_bonds
 304                                 #       atoms.each { |ap|
 305                                 #               j = ap.index
 306                                 #               (j + 1 ... natoms).each { |k|
 307                                 #                       if calc_bond(j, k) < 1.7
 308                                 #                               create_bond(j, k)
 309                                 #                       end
 310                                 #               }
 311                                 #       }
 312                                         new_unit = false
 313                                 #       create_frame
 314                                 else
 315                                     dont_create = false
 316                                         if (search_mode == 1 && nsearch == 1) || first_line
 317                                                 #  The input coordinate and the first frame for geometry search
 318                                                 #  can have the same coordinate as the last frame; if this is the case, then
 319                                                 #  do not create the new frame
 320                                                 select_frame(nframes - 1)
 321                                                 dont_create = true
 322                                                 each_atom { |ap|
 323                                                   if (ap.r - coords[ap.index]).length2 > 1e-8
 324                                                     dont_create = false
 325                                                         break
 326                                                   end
 327                                                 }
 328                                         end
 329                                         if !dont_create
 330                                                 create_frame([coords])  #  Should not be (coords)
 331                                         end
 332                                 end
 333                                 set_property("energy", energy) if energy
 334                         elsif line =~ /BEGINNING GEOMETRY SEARCH POINT/
 335                                 energy = nil   #  New search has begun, so clear the energy
 336                                 search_mode = 1
 337                         elsif line =~ /CONSTRAINED OPTIMIZATION POINT/
 338                                 energy = nil   #  New search has begun, so clear the energy
 339                                 search_mode = 2
 340                         elsif line =~ /FINAL .* ENERGY IS *([-.0-9]+) AFTER/
 341                                 if search_mode != 2
 342                                         energy = $1.to_f
 343                                         set_property("energy", energy)
 344                                 end
 345                         elsif line =~ /TOTAL ENERGY += +([-.0-9]+)/
 346                                 energy = $1.to_f
 347                         elsif false && line =~ /EQUILIBRIUM GEOMETRY LOCATED/i
 348                                 set_progress_message(mes + "\nReading optimized coordinates...")
 349                                 fp.gets; fp.gets; fp.gets
 350                                 n = 0
 351                                 while (line = fp.gets) != nil
 352                                         break if line =~ /^\s*$/
 353                                         line.chomp
 354                                         atom, an, x, y, z = line.split
 355                                         ap = atoms[n]
 356                                         ap.r = Vector3D[x, y, z]
 357                                         n += 1
 358                                         break if n >= natoms
 359                                 end
 360                         #       if ne_alpha > 0 && ne_beta > 0
 361                         #               #  Allocate basis set record again, to update the atomic coordinates
 362                         #               allocate_basis_set_record(rflag, ne_alpha, ne_beta)
 363                         #       end
 364                         elsif line =~ /ATOMIC BASIS SET/
 365                                 lines = []
 366                                 lineno = fp.lineno
 367                                 while (line = fp.gets)
 368                                         break if line =~ /TOTAL NUMBER OF BASIS SET/
 369                                         line.chomp!
 370                                         lines.push(line)
 371                                 end
 372                                 ncomps = sub_load_gamess_log_basis_set(lines, lineno)
 373                         elsif line =~ /NUMBER OF OCCUPIED ORBITALS/
 374                                 line =~ /=\s*(\d+)/
 375                                 n = Integer($1)
 376                                 if line =~ /ALPHA/
 377                                         ne_alpha = n
 378                                 else
 379                                         ne_beta = n
 380                                 end
 381                         elsif line =~ /SCFTYP=(\w+)/
 382                                 scftyp = $1
 383                                 if ne_alpha > 0 || ne_beta > 0
 384                                         rflag = 0
 385                                         case scftyp
 386                                         when "RHF"
 387                                                 rflag = 1
 388                                         when "ROHF"
 389                                                 rflag = 2
 390                                         end
 391                                 end
 392                         elsif line =~ /(ALPHA|BETA)\s*SET/
 393                                 alpha_beta = $1
 394                         elsif line =~ /^\s*(EIGENVECTORS|MOLECULAR ORBITALS)\s*$/
 395                                 if mo_count == 0
 396                                         clear_mo_coefficients
 397                                         set_mo_info(:type=>["UHF", "RHF", "ROHF"][rflag], :alpha=>ne_alpha, :beta=>ne_beta)
 398                                 end
 399                                 mo_count += 1
 400                                 line = fp.gets; line = fp.gets
 401                                 lineno = fp.lineno
 402                                 lines = []
 403                                 set_progress_message(mes + "\nReading MO coefficients...")
 404                                 while (line = fp.gets)
 405                                         break if line =~ /\.\.\.\.\.\./ || line =~ /----------------/
 406                                         line.chomp!
 407                                         lines.push(line)
 408                                 end
 409                                 sub_load_gamess_log_mo_coefficients(lines, lineno, ncomps)
 410                                 set_progress_message(mes)
 411                         elsif line =~ /N A T U R A L   B O N D   O R B I T A L   A N A L Y S I S/
 412                                 nbo_lines = []
 413                                 while (line = fp.gets) != nil
 414                                   break if line =~ /done with NBO analysis/
 415                                   nbo_lines.push(line)
 416                                 end
 417                                 import_nbo_log(nbo_lines)
 418                         end
 419                 end
 420         end
 421         if nframes > 0
 422           select_frame(nframes - 1)
 423         end
 424         if energy && energy != 0.0
 425           set_property("energy", energy)
 426         end
 427         hide_progress_panel
 428         self.update_enabled = true
 429         (n > 0 ? true : false)
 430   end
 431
 432   def sub_load_gaussian_log(fp)
 433
 434     if natoms == 0
 435                 new_unit = true
 436         else
 437                 new_unit = false
 438     end
 439         if nframes > 0
 440                 create_frame
 441                 frame = nframes - 1
 442         end
 443         n = 0
 444         nf = 0
 445         energy = nil
 446         use_input_orientation = false
 447         show_progress_panel("Loading Gaussian out file...")
 448         while 1
 449                 line = fp.gets
 450                 if line == nil
 451                         break
 452                 end
 453                 line.chomp!
 454                 if line =~ /(Input|Standard) orientation/
 455                         match = $1
 456                         if match == "Input"
 457                                 use_input_orientation = true if nf == 0
 458                                 next if !use_input_orientation
 459                         else
 460                                 next if use_input_orientation
 461                         end
 462                         4.times { line = fp.gets }    #  Skip four lines
 463                         n = 0
 464                         coords = []
 465                         anums = []
 466                         while (line = fp.gets) != nil
 467                                 break if line =~ /-----/
 468                                 num, charge, type, x, y, z = line.split
 469                                 coords.push(Vector3D[x, y, z])
 470                                 anums.push(charge)
 471                                 n += 1
 472                         end
 473                         if new_unit
 474                                 #  Build a new molecule
 475                                 anums.each_index { |i|
 476                                         ap = add_atom("X")
 477                                         ap.atomic_number = anums[i]
 478                                         ap.atom_type = ap.element
 479                                         ap.name = sprintf("%s%d", ap.element, i)
 480                                         ap.r = coords[i]
 481                                 }
 482                                 #  Find bonds
 483                         #       atoms.each { |ap|
 484                         #               j = ap.index
 485                         #               (j + 1 ... natoms).each { |k|
 486                         #                       if calc_bond(j, k) < 1.7
 487                         #                               create_bond(j, k)
 488                         #                       end
 489                         #               }
 490                         #       }
 491                                 guess_bonds
 492                                 new_unit = false
 493                         #       create_frame
 494                         else
 495                                 create_frame([coords])  #  Should not be (coords)
 496                         end
 497                         if energy
 498                                 # TODO: to ensure whether the energy output line comes before
 499                                 # or after the atomic coordinates.
 500                                 set_property("energy", energy)
 501                         end
 502                         nf += 1
 503                 elsif line =~ /SCF Done: *E\(\w+\) *= *([-.0-9]+)/
 504                         energy = $1.to_f
 505                 end
 506         end
 507         if energy
 508                 set_property("energy", energy)
 509         end
 510         hide_progress_panel
 511         (n > 0 ? true : false)
 512   end
 513
 514   def sub_load_psi4_log(fp)
 515     if natoms == 0
 516       new_unit = true
 517     else
 518       new_unit = false
 519     end
 520     n = 0
 521     nf = 0
 522     energy = nil
 523
 524     show_progress_panel("Loading Psi4 output file...")
 525
 526     getline = lambda { @lineno += 1; return fp.gets }
 527
 528     #  Import coordinates and energies
 529     vecs = []
 530     ats = []
 531     first_frame = nframes
 532     trans = nil
 533     hf_type = nil
 534     nalpha = nil
 535     nbeta = nil
 536     while line = getline.call
 537       if line =~ /==> Geometry <==/
 538         #  Skip until line containing "------"
 539         while line = getline.call
 540           break if line =~ /------/
 541         end
 542         vecs.clear
 543         index = 0
 544         #  Read atom positions
 545         while line = getline.call
 546           line.chomp!
 547           break if line =~ /^\s*$/
 548           tokens = line.split(' ')
 549           if natoms > 0 && first_frame == nframes
 550             if index >= natoms || tokens[0] != atoms[index].element
 551               hide_progress_panel
 552               raise MolbyError, "The atom list does not match the current structure at line #{@lineno}"
 553             end
 554           end
 555           vecs.push(Vector3D[Float(tokens[1]), Float(tokens[2]), Float(tokens[3])])
 556           if natoms == 0
 557             ats.push(tokens[0])
 558           end
 559           index += 1
 560         end
 561         if natoms == 0
 562           #  Create molecule from the initial geometry
 563           ats.each_with_index { |aname, i|
 564             #  Create atoms
 565             ap = add_atom(aname)
 566             ap.element = aname
 567             ap.atom_type = ap.element
 568             ap.name = sprintf("%s%d", aname, i)
 569             ap.r = vecs[i]
 570           }
 571           guess_bonds
 572         else
 573           if vecs.length != natoms
 574             break  #  Log file is incomplete
 575           end
 576           #  Does this geometry differ from the last one?
 577           vecs.length.times { |i|
 578             if (atoms[i].r - vecs[i]).length2 > 1.0e-14
 579               #  Create a new frame and break
 580               create_frame
 581               vecs.length.times { |j|
 582                 atoms[j].r = vecs[j]
 583               }
 584               break
 585             end
 586           }
 587         end
 588         #  end geometry
 589       elsif line =~ /Final Energy: +([-.0-9]+)/
 590         #  Energy for this geometry
 591         energy = Float($1)
 592         set_property("energy", energy)
 593         if line =~ /RHF/
 594           hf_type = "RHF"
 595         elsif line =~ /UHF/
 596           hf_type = "UHF"
 597         elsif line =~ /ROHF/
 598           hf_type = "ROHF"
 599         end
 600       elsif line =~ /^ *Nalpha *= *(\d+)/
 601         nalpha = Integer($1)
 602       elsif line =~ /^ *Nbeta *= *(\d+)/
 603         nbeta = Integer($1)
 604       end
 605     end
 606     hide_progress_panel
 607     clear_basis_set
 608     clear_mo_coefficients
 609     set_mo_info(:type => hf_type, :alpha => nalpha, :beta => nbeta)
 610     return true
 611   end
 612
 613   #  mol.set_mo_info should be set before calling this function
 614   def sub_load_molden(fp)
 615     getline = lambda { @lineno += 1; return fp.gets }
 616     bohr = 0.529177210903
 617     errmsg = nil
 618     ncomps = 0  #  Number of components (AOs)
 619     occ_alpha = 0  #  Number of occupied alpha orbitals
 620     occ_beta = 0   #  Number of occupied beta orbitals
 621     catch :ignore do
 622       while line = getline.call
 623         if line =~ /^\[Atoms\]/
 624           i = 0
 625           while line = getline.call
 626             if line =~ /^[A-Z]/
 627               #  element, index, atomic_number, x, y, z (in AU)
 628               a = line.split(' ')
 629               if atoms[i].atomic_number != Integer(a[2]) ||
 630                 (atoms[i].x - Float(a[3]) * bohr).abs > 1e-4 ||
 631                 (atoms[i].y - Float(a[4]) * bohr).abs > 1e-4 ||
 632                 (atoms[i].z - Float(a[5]) * bohr).abs > 1e-4
 633                 errmsg = "The atom list does not match the current molecule."
 634                 throw :ignore
 635               end
 636               i += 1
 637             else
 638               break
 639             end
 640           end
 641           redo  #  The next line will be the beginning of the next block
 642         elsif line =~ /^\[GTO\]/
 643           shell = 0
 644           atom_index = 0
 645           while line = getline.call
 646             #  index, 0?
 647             a = line.split(' ')
 648             break if a.length != 2
 649             #  loop for shells
 650             while line = getline.call
 651               #  type, no_of_primitives, 1.00?
 652               a = line.split(' ')
 653               break if a.length != 3   #  Terminated by a blank line
 654               case a[0]
 655               when "s"
 656                 sym = 0; n = 1
 657               when "p"
 658                 sym = 1; n = 3
 659               when "d"
 660                 sym = 2; n = 6    #  TODO: handle both spherical and cartesian
 661               when "f"
 662                 sym = 3; n = 10
 663               when "g"
 664                 sym = 4; n = 15
 665               else
 666                 raise MolbyError, "Unknown gaussian shell type '#{a[0]}' at line #{@lineno} in MOLDEN file"
 667               end
 668               nprimitives = Integer(a[1])
 669               nprimitives.times { |i|
 670                 line = getline.call   #  exponent, contraction
 671                 b = line.split(' ')
 672                 add_gaussian_primitive_coefficients(Float(b[0]), Float(b[1]), 0.0)
 673               }
 674               #  end of one shell
 675               add_gaussian_orbital_shell(atom_index, sym, nprimitives)
 676               shell += 1
 677               ncomps += n
 678             end
 679             #  end of one atom
 680             atom_index += 1
 681           end
 682           redo  #  The next line will be the beginning of the next block
 683         elsif line =~ /^\[MO\]/
 684           m = []
 685           idx_alpha = 1   #  set_mo_coefficients() accepts 1-based index of MO
 686           idx_beta = 1
 687           while true
 688             #  Loop for each MO
 689             m.clear
 690             ene = nil
 691             spin = nil
 692             sym = nil   #  Not used in Molby
 693             occ = nil
 694             i = 0
 695             while line = getline.call
 696               if line =~ /^ *Sym= *(\w+)/
 697                 sym = $1
 698               elsif line =~ /^ *Ene= *([-+.0-9e]+)/
 699                 ene = Float($1)
 700               elsif line =~ /^ *Spin= *(\w+)/
 701                 spin = $1
 702               elsif line =~ /^ *Occup= *([-+.0-9e]+)/
 703                 occ = Float($1)
 704                 if occ > 0.0
 705                   if spin == "Alpha"
 706                     occ_alpha += 1
 707                   else
 708                     occ_beta += 1
 709                   end
 710                 end
 711               elsif line =~ /^ *([0-9]+) +([-+.0-9e]+)/
 712                 m[i] = Float($2)
 713                 i += 1
 714                 if i >= ncomps
 715                   if spin == "Alpha"
 716                     idx = idx_alpha
 717                     idx_alpha += 1
 718                   else
 719                     idx = idx_beta
 720                     idx_beta += 1
 721                   end
 722                   set_mo_coefficients(idx, ene, m)
 723                   break
 724                 end
 725               else
 726                 break
 727               end
 728             end
 729             break if i < ncomps  #  no MO info was found
 730           end
 731           next
 732         end
 733       end   #  end while
 734     end     #  end catch
 735     if errmsg
 736       message_box("The MOLDEN file was found but not imported. " + errmsg, "Psi4 import info", :ok)
 737       return false
 738     end
 739     return true
 740   end
 741
 742   def loadout(filename)
 743   retval = false
 744   fp = open(filename, "rb")
 745   @lineno = 0
 746   begin
 747     while s = fp.gets
 748       @lineno += 1
 749       if s =~ /Gaussian/
 750         retval = sub_load_gaussian_log(fp)
 751         break
 752       elsif s =~ /GAMESS/
 753         retval = sub_load_gamess_log(fp)
 754         break
 755       elsif s =~ /Psi4/
 756         retval = sub_load_psi4_log(fp)
 757         if retval
 758           #  If .molden file exists, then try to read it
 759           mname = filename.gsub(/\.\w*$/, ".molden")
 760           if File.exists?(mname)
 761             fp2 = open(mname, "rb")
 762             if fp2
 763               flag = sub_load_molden(fp2)
 764               fp2.close
 765             end
 766           end
 767         end
 768         break
 769       end
 770     end
 771   rescue
 772     hide_progress_panel
 773     raise
 774   end
 775         fp.close
 776         return retval
 777   end
 778
 779   alias :loadlog :loadout
 780
 781   def loadxyz(filename)
 782         fp = open(filename, "rb")
 783         n = 0
 784         coords = []
 785         names = []
 786         cell = nil
 787         while 1
 788           line = fp.gets
 789           if line == nil
 790                 fp.close
 791                 break
 792           end
 793           line.chomp
 794           toks = line.split
 795           if coords.length == 0
 796             #  Allow "number of atoms" line or "number of atoms and crystallographic cell parameter"
 797                 #  (Chem3D xyz format)
 798                 next if toks.length == 1
 799                 if toks.length == 7
 800                   cell = toks[1..6].map { |s| Float(s.sub(/\(\d+\)/, "")) }  #  Remove (xx) and convert to float
 801                   next
 802                 end
 803           end
 804           name, x, y, z = line.split
 805           next if z == nil
 806           x = Float(x.sub(/\(\d+\)/, ""))
 807           y = Float(y.sub(/\(\d+\)/, ""))
 808           z = Float(z.sub(/\(\d+\)/, ""))
 809           r = Vector3D[x, y, z]
 810           coords.push(r)
 811           names.push(name)
 812           n += 1
 813         end
 814         celltr = nil
 815         if cell
 816           self.cell = cell
 817           celltr = self.cell_transform
 818         end
 819         names.each_index { |i|
 820           ap = add_atom(names[i])
 821           names[i] =~ /^([A-Za-z]{1,2})/
 822           element = $1.capitalize
 823           ap.element = element
 824           ap.atom_type = element
 825           if celltr
 826             ap.r = celltr * coords[i]
 827           else
 828             ap.r = coords[i]
 829           end
 830         }
 831         guess_bonds
 832         #  Find bonds
 833 #       atoms.each { |ap|
 834 #         j = ap.index
 835 #         (j + 1 ... natoms).each { |k|
 836 #               if calc_bond(j, k) < 1.7
 837 #               #  create_bond(j, k)
 838 #               end
 839 #         }
 840 #       }
 841 #       self.undo_enabled = save_undo_enabled
 842         (n > 0 ? true : false)
 843   end
 844
 845   def savexyz(filename)
 846     open(filename, "wb") { |fp|
 847           fp.printf "%d\n", self.natoms
 848           each_atom { |ap|
 849             fp.printf "%s %.5f %.5f %.5f\n", ap.element, ap.x, ap.y, ap.z
 850           }
 851         }
 852         return true
 853   end
 854
 855   def loadzmat(filename)
 856     self.remove(All)
 857         open(filename, "rb") { |fp|
 858           while (line = fp.gets)
 859             line.chomp!
 860                 a = line.split
 861                 an = Molecule.guess_atomic_number_from_name(a[0])
 862                 elm = Parameter.builtin.elements[an].name
 863                 base1 = a[1].to_i - 1
 864                 base2 = a[3].to_i - 1
 865                 base3 = a[5].to_i - 1
 866                 base1 = nil if base1 < 0
 867                 base2 = nil if base2 < 0
 868                 base3 = nil if base3 < 0
 869                 add_atom(a[0], elm, elm, a[2].to_f, base1, a[4].to_f, base2, a[6].to_f, base3)
 870           end
 871         }
 872         return true
 873   end
 874
 875   def savezmat(filename)
 876   end
 877
 878   def loadcom(filename)
 879 #       save_undo_enabled = self.undo_enabled?
 880 #       self.undo_enabled = false
 881         self.remove(All)
 882         fp = open(filename, "rb")
 883         section = 0
 884         while (line = fp.gets)
 885           line.chomp!
 886           if section == 0
 887             section = 1 if line =~ /^\#/
 888                 next
 889           elsif section == 1 || section == 2
 890             section += 1 if line =~ /^\s*$/
 891                 next
 892           else
 893             #  The first line is skipped (charge and multiplicity)
 894                 while (line = fp.gets)
 895                   line.chomp!
 896                   break if line =~ /^\s*$/
 897                   a = line.split(/\s*[ \t,\/]\s*/)
 898                   r = Vector3D[Float(a[1]), Float(a[2]), Float(a[3])]
 899                   ap = add_atom(a[0])
 900                   a[0] =~ /^([A-Za-z]{1,2})/
 901                   element = $1.capitalize
 902               ap.element = element
 903               ap.atom_type = element
 904               ap.r = r
 905                 end
 906                 break
 907           end
 908         end
 909         fp.close
 910         guess_bonds
 911 #       self.undo_enabled = save_undo_enabled
 912         return true
 913   end
 914
 915   def loadinp(filename)
 916 #       save_undo_enabled = self.undo_enabled?
 917 #       self.undo_enabled = false
 918         self.remove(All)
 919         fp = open(filename, "rb")
 920         section = 0
 921         has_basis = false
 922         while (line = fp.gets)
 923           if line =~ /\A \$BASIS/
 924             has_basis = true
 925                 next
 926           end
 927           next if line !~ /\A \$DATA/
 928           line = fp.gets   #  Title line
 929           line = fp.gets   #  Symmetry line
 930           while (line = fp.gets)
 931 #           puts line
 932             line.chomp!
 933                 break if line =~ /\$END/
 934                 a = line.split
 935                 ap = add_atom(a[0])
 936                 ap.atomic_number = Integer(a[1])
 937                 ap.atom_type = ap.element
 938                 r = Vector3D[Float(a[2]), Float(a[3]), Float(a[4])]
 939                 ap.r = r;
 940                 if !has_basis
 941                   #  Skip until one blank line is detected
 942                   while (line = fp.gets) && line =~ /\S/
 943                   end
 944                 end
 945       end
 946           break
 947         end
 948         fp.close
 949         guess_bonds
 950 #       self.undo_enabled = save_undo_enabled
 951         return true
 952   end
 953
 954   def saveinp(filename)
 955     if natoms == 0
 956       raise MolbyError, "cannot save GAMESS input; the molecule is empty"
 957     end
 958     fp = open(filename, "wb")
 959         now = Time.now.to_s
 960         fp.print <<end_of_header
 961 !  GAMESS input
 962 !  Generated by Molby at #{now}
 963  $CONTRL COORD=UNIQUE EXETYP=RUN ICHARG=0
 964          ICUT=20 INTTYP=HONDO ITOL=30
 965          MAXIT=200 MOLPLT=.T. MPLEVL=0
 966          MULT=1 QMTTOL=1e-08 RUNTYP=OPTIMIZE
 967          SCFTYP=RHF UNITS=ANGS                  $END
 968  $SCF    CONV=1.0E-06 DIRSCF=.T.                $END
 969  $STATPT NSTEP=400 OPTTOL=1.0E-06               $END
 970  $SYSTEM MEMDDI=0 MWORDS=16 TIMLIM=50000        $END
 971  $BASIS  GBASIS=N31 NDFUNC=1 NGAUSS=6           $END
 972  $GUESS  GUESS=HUCKEL                           $END
 973 !
 974  $DATA
 975  #{name}
 976  C1
 977 end_of_header
 978         each_atom { |ap|
 979                 next if ap.atomic_number == 0
 980                 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
 981         }
 982         fp.print " $END\n"
 983         fp.close
 984         return true
 985   end
 986
 987   def savecom(filename)
 988     if natoms == 0
 989       raise MolbyError, "cannot save Gaussian input; the molecule is empty"
 990     end
 991     fp = open(filename, "wb")
 992         base = File.basename(filename, ".*")
 993         fp.print <<end_of_header
 994 %Chk=#{base}.chk
 995 \# PM3 Opt
 996
 997  #{name}; created by Molby at #{Time.now.to_s}
 998
 999  0 1
1000 end_of_header
1001         each_atom { |ap|
1002             next if ap.atomic_number == 0
1003                 fp.printf "%-6s %10.6f %10.6f %10.6f\n", ap.element, ap.r.x, ap.r.y, ap.r.z
1004         }
1005         fp.print "\n"
1006         fp.close
1007         return true
1008   end
1009
1010   alias :loadgjf :loadcom
1011   alias :savegjf :savecom
1012
1013   def loadcif(filename)
1014     mol = self
1015     def getciftoken(fp)
1016           while @tokens.length == 0
1017             line = fp.gets
1018             return nil if !line
1019                 if line[0] == ?;
1020                   s = line
1021                   while line = fp.gets
1022                     break if line[0] == ?;
1023                     s += line
1024                   end
1025                   return s if !line
1026                   s += ";"
1027                   @tokens.push(s)
1028                   line = line[1...line.length]
1029                 end
1030             line.strip!
1031             line.scan(/'[^\']*'|"[^\"]*"|[^#\s]+|#.*/) do |s|
1032               next if s[0] == ?#
1033                   if s =~ /^data_|loop_|global_|save_|stop_/i
1034                     s = "#" + s    #  Label for reserved words
1035                   end
1036                   @tokens.push(s)
1037             end
1038           end
1039           return @tokens.shift
1040         end
1041         def float_strip_rms(str)
1042           str =~ /^(-?)(\d*)(\.(\d*))?(\((\d+)\))?$/
1043           sgn, i, frac, rms = $1, $2, $4, $6
1044           i = i.to_f
1045           if frac
1046                 base = 0.1 ** frac.length
1047                 i = i + frac.to_f * base
1048           else
1049             base = 1.0
1050           end
1051           if rms
1052             rms = rms.to_f * base
1053           else
1054             rms = 0.0
1055           end
1056           if sgn == "-"
1057             i = -i
1058           end
1059           return i, rms
1060         end
1061         def parse_symmetry_operation(str)
1062           if str == "."
1063             sym = nil
1064           elsif (str =~ /(\d+)_(\d)(\d)(\d)/) || (str =~ /(\d+) +(\d)(\d)(\d)/)
1065             sym = [Integer($1) - 1, Integer($2) - 5, Integer($3) - 5, Integer($4) - 5]
1066           elsif (str =~ /^(\d+)$/)
1067             sym = [Integer($1) - 1, 0, 0, 0]
1068           end
1069           if sym && (sym[0] == 0 && sym[1] == 0 && sym[2] == 0 && sym[3] == 0)
1070             sym = nil
1071           end
1072           sym
1073         end
1074         def find_atom_by_name(mol, name)
1075           name = name.delete(" ()")
1076           ap = mol.atoms[name] rescue ap = nil
1077           return ap
1078         end
1079         selfname = self.name
1080         fp = open(filename, "rb")
1081         data_identifier = nil
1082         @tokens = []
1083         count_up = 1
1084         while true
1085           warn_message = ""
1086           verbose = nil
1087           bond_defined = false
1088           special_positions = []
1089           mol.remove(All)
1090           cell = []
1091           cell_trans = cell_trans_inv = Transform.identity
1092           token = getciftoken(fp)
1093           pardigits_re = /\(\d+\)/
1094           calculated_atoms = []
1095           while token != nil
1096             if token =~ /^\#data_/i
1097                   if data_identifier == nil || mol.natoms == 0
1098                     #  First block or no atoms yet
1099             #  Continue processing of this molecule
1100                     data_identifier = token
1101                         token = getciftoken(fp)
1102                         next
1103                   else
1104                     #  Description of another molecule begins here
1105                         data_identifier = token
1106                         break
1107                   end
1108             elsif token =~ /^_cell/
1109                   val = getciftoken(fp)
1110                   if token == "_cell_length_a"
1111                     cell[0], cell[6] = float_strip_rms(val)
1112                   elsif token == "_cell_length_b"
1113                     cell[1], cell[7] = float_strip_rms(val)
1114                   elsif token == "_cell_length_c"
1115                     cell[2], cell[8] = float_strip_rms(val)
1116                   elsif token == "_cell_angle_alpha"
1117                     cell[3], cell[9] = float_strip_rms(val)
1118                   elsif token == "_cell_angle_beta"
1119                     cell[4], cell[10] = float_strip_rms(val)
1120                   elsif token == "_cell_angle_gamma"
1121                     cell[5], cell[11] = float_strip_rms(val)
1122                   end
1123                   if cell.length == 12 && cell.all?
1124                     mol.cell = cell
1125             puts "Unit cell is set to #{mol.cell.inspect}." if verbose
1126                     cell = []
1127                     cell_trans = mol.cell_transform
1128                     cell_trans_inv = cell_trans.inverse
1129                   end
1130                   token = getciftoken(fp)
1131                   next
1132         elsif token.casecmp("#loop_") == 0
1133               labels = []
1134                   while (token = getciftoken(fp)) && token[0] == ?_
1135                     labels.push(token)
1136                   end
1137                   if labels[0] =~ /symmetry_equiv_pos|space_group_symop|atom_site_label|atom_site_aniso_label|geom_bond/
1138                     hlabel = Hash.new(-10000000)
1139                     labels.each_with_index { |lb, i|
1140                           hlabel[lb] = i
1141                     }
1142                     data = []
1143                     n = labels.length
1144                     a = []
1145                     while true
1146                           break if token == nil || token[0] == ?_ || token[0] == ?#
1147                           a.push(token)
1148                           if a.length == n
1149                             data.push(a)
1150                             a = []
1151                           end
1152                           token = getciftoken(fp)
1153                     end
1154                     if labels[0] =~ /^_symmetry_equiv_pos/ || labels[0] =~ /^_space_group_symop/
1155                       data.each { |d|
1156                             symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]] || d[hlabel["_space_group_symop_operation_xyz"]]
1157                             symstr.delete("\"\'")
1158                             exps = symstr.split(/,/)
1159                             sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
1160                             exps.each_with_index { |s, i|
1161                               terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
1162                                   terms.each { |a|
1163                                     #  a[0]: sign, a[2]: numerator, a[4]: denometer
1164                                     if a[4] != nil
1165                                       #  The number part is a[2]/a[4]
1166                                       num = Float(a[2])/Float(a[4])
1167                                     elsif a[2] != nil
1168                                       #  The number part is either integer or a floating point
1169                                       num = Float(a[2])
1170                                     else
1171                                       num = 1.0
1172                                     end
1173                                     num = -num if a[0][0] == ?-
1174                                     xyz = (a[5] || a[6])
1175                                     if xyz == "x" || xyz == "X"
1176                                       sym[i] = num
1177                                     elsif xyz == "y" || xyz == "Y"
1178                                       sym[i + 3] = num
1179                                     elsif xyz == "z" || xyz == "Z"
1180                                       sym[i + 6] = num
1181                                     else
1182                                       sym[9 + i] = num
1183                                     end
1184                                   }
1185                             }
1186                             puts "symmetry operation #{sym.inspect}" if verbose
1187                             mol.add_symmetry(Transform.new(sym))
1188                           }
1189                           puts "#{mol.nsymmetries} symmetry operations are added" if verbose
1190                     elsif labels[0] =~ /^_atom_site_label/
1191                           #  Create atoms
1192                           data.each { |d|
1193                             name = d[hlabel["_atom_site_label"]]
1194                             elem = d[hlabel["_atom_site_type_symbol"]]
1195                             fx = d[hlabel["_atom_site_fract_x"]]
1196                             fy = d[hlabel["_atom_site_fract_y"]]
1197                             fz = d[hlabel["_atom_site_fract_z"]]
1198                             uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
1199                             biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
1200                             occ = d[hlabel["_atom_site_occupancy"]]
1201                             calc = d[hlabel["_atom_site_calc_flag"]]
1202                             name = name.delete(" ()")
1203                             if elem == nil || elem == ""
1204                               if name =~ /[A-Za-z]{1,2}/
1205                                     elem = $&.capitalize
1206                                   else
1207                                     elem = "Du"
1208                                   end
1209                             end
1210                             ap = mol.add_atom(name, elem, elem)
1211                             ap.fract_x, ap.sigma_x = float_strip_rms(fx)
1212                             ap.fract_y, ap.sigma_y = float_strip_rms(fy)
1213                             ap.fract_z, ap.sigma_z = float_strip_rms(fz)
1214                             if biso
1215                               ap.temp_factor, sig = float_strip_rms(biso)
1216                             elsif uiso
1217                               ap.temp_factor, sig = float_strip_rms(uiso)
1218                                   ap.temp_factor *= 78.9568352087149          #  8*pi*pi
1219                             end
1220                             ap.occupancy, sig = float_strip_rms(occ)
1221                             if calc == "c" || calc == "calc"
1222                               calculated_atoms.push(ap.index)
1223                         end
1224                             #  Guess special positions
1225                             (1...mol.nsymmetries).each { |isym|
1226                               sr = ap.fract_r
1227                               sr = (mol.transform_for_symop(isym) * sr) - sr;
1228                                   nx = (sr.x + 0.5).floor
1229                                   ny = (sr.y + 0.5).floor
1230                                   nz = (sr.z + 0.5).floor
1231                                   if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
1232                                     #  [isym, -nx, -ny, -nz] transforms this atom to itself
1233                                     #  The following line is equivalent to:
1234                                     #    if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
1235                                     #    special_positions[ap.index].push(...)
1236                                     (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
1237                                   end
1238                             }
1239                             if verbose && special_positions[ap.index]
1240                               puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
1241                             end
1242                           }
1243                           puts "#{mol.natoms} atoms are created." if verbose
1244                     elsif labels[0] =~ /^_atom_site_aniso_label/
1245                       #  Set anisotropic parameters
1246                           c = 0
1247                           data.each { |d|
1248                             name = d[hlabel["_atom_site_aniso_label"]]
1249                             ap = find_atom_by_name(mol, name)
1250                             next if !ap
1251                             u11 = d[hlabel["_atom_site_aniso_U_11"]]
1252                             if u11
1253                               usig = []
1254                               u11, usig[0] = float_strip_rms(u11)
1255                               u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
1256                               u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
1257                               u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
1258                               u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
1259                               u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
1260                               ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
1261                                   c += 1
1262                             end
1263                           }
1264                           puts "#{c} anisotropic parameters are set." if verbose
1265                     elsif labels[0] =~ /^_geom_bond/
1266                       #  Create bonds
1267                           exbonds = []
1268                           data.each { |d|
1269                             n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
1270                             n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
1271                             sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
1272                             sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
1273                             n1 = find_atom_by_name(mol, n1)
1274                             n2 = find_atom_by_name(mol, n2)
1275                             next if n1 == nil || n2 == nil
1276                         n1 = n1.index
1277                             n2 = n2.index
1278                             sym1 = parse_symmetry_operation(sym1)
1279                             sym2 = parse_symmetry_operation(sym2)
1280                             if sym1 || sym2
1281                               exbonds.push([n1, n2, sym1, sym2])
1282                             else
1283                               mol.create_bond(n1, n2)
1284                             end
1285                             tr1 = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
1286                             tr2 = (sym2 ? mol.transform_for_symop(sym2) : Transform.identity)
1287                             if special_positions[n1]
1288                                   #  Add extra bonds for equivalent positions of n1
1289                                   special_positions[n1].each { |symop|
1290                                     sym2x = mol.symop_for_transform(tr1 * mol.transform_for_symop(symop) * tr1.inverse * tr2)
1291                                     exbonds.push([n1, n2, sym1, sym2x])
1292                                   }
1293                             end
1294                             if special_positions[n2]
1295                                   #  Add extra bonds n2-n1.symop, where symop transforms n2 to self
1296                                   tr = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
1297                                   special_positions[n2].each { |symop|
1298                                     sym1x = mol.symop_for_transform(tr2 * mol.transform_for_symop(symop) * tr2.inverse * tr1)
1299                                     exbonds.push([n2, n1, sym2, sym1x])
1300                                   }
1301                             end
1302                       }
1303               if mol.nbonds > 0
1304                 bond_defined = true
1305               end
1306                           puts "#{mol.nbonds} bonds are created." if verbose
1307                           if calculated_atoms.length > 0
1308                             #  Guess bonds for calculated hydrogen atoms
1309                             n1 = 0
1310                             calculated_atoms.each { |ai|
1311                               if mol.atoms[ai].connects.length == 0
1312                                     as = mol.find_close_atoms(ai)
1313                                     as.each { |aj|
1314                                       mol.create_bond(ai, aj)
1315                                           n1 += 1
1316                                     }
1317                                   end
1318                             }
1319                             puts "#{n1} bonds are guessed." if verbose
1320                           end
1321                           if exbonds.length > 0
1322                             h = Dialog.run("CIF Import: Symmetry Expansion") {
1323                               layout(1,
1324                                     item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
1325                                     item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
1326                                     item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
1327                                     item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
1328                                   )
1329                                   radio_group("atoms_only", "fragment", "ignore")
1330                             }
1331                             if h[:status] == 0 && h["ignore"] == 0
1332                               atoms_only = (h["atoms_only"] != 0)
1333                                   if !atoms_only
1334                                     fragments = []
1335                                     mol.each_fragment { |f| fragments.push(f) }
1336                                   end
1337                                   debug = nil
1338                                   exbonds.each { |ex|
1339                                     if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
1340                                     ex0 = ex.dup
1341                                     (2..3).each { |i|
1342                                       symop = ex[i]
1343                                           if symop == nil
1344                                             ex[i + 2] = ex[i - 2]
1345                                           else
1346                                             if debug; puts "  symop = #{symop.inspect}"; end
1347                                             #  Expand the atom or the fragment including the atom
1348                                             if atoms_only
1349                                                   ig = IntGroup[ex[i - 2]]
1350                                                   idx = 0
1351                                             else
1352                                                   ig = fragments.find { |f| f.include?(ex[i - 2]) }
1353                                                   ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
1354                                             end
1355                                             symop[4] = ex[i - 2]  #  Base atom
1356                                             if debug; puts "  expanding #{ig} by #{symop.inspect}"; end
1357                                             a = mol.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
1358                                             ex[i + 2] = a[idx]   #  Index of the expanded atom
1359                                           end
1360                                     }
1361                                     if ex[4] && ex[5] && ex[4] != ex[5]
1362                                       if debug; puts "  creating bond #{ex[4]} - #{ex[5]}"; end
1363                                       mol.create_bond(ex[4], ex[5])
1364                                     end
1365                                   }
1366                             end
1367                           end
1368                           puts "#{mol.nbonds} bonds are created." if verbose
1369                     end
1370                     next
1371                   else
1372                   #  puts "Loop beginning with #{labels[0]} is skipped"
1373                   end
1374             else
1375               #  Skip this token
1376                   token = getciftoken(fp)
1377             end
1378             #  Skip tokens until next tag or reserved word is detected
1379             while token != nil && token[0] != ?_ && token[0] != ?#
1380                   token = getciftoken(fp)
1381             end
1382             next
1383           end
1384       if !bond_defined
1385                 mol.guess_bonds
1386           end
1387           if token != nil && token == data_identifier
1388             #  Process next molecule: open a new molecule and start adding atom on that
1389                 mol = Molecule.new
1390                 count_up += 1
1391                 (@aux_mols ||= []).push(mol)
1392                 next
1393           end
1394           break
1395         end
1396         fp.close
1397 #       self.undo_enabled = save_undo_enabled
1398         return true
1399   end
1400
1401   def dump(group = nil)
1402     def quote(str)
1403           if str == ""
1404             str = "\"\""
1405           else
1406             str = str.gsub(/%/, "%%")
1407                 str.gsub!(/ /, "%20")
1408                 str.gsub!(/\t/, "%09")
1409           end
1410           str
1411         end
1412         group = atom_group(group ? group : 0...natoms)
1413         s = ""
1414         group.each { |i|
1415           ap = atoms[i]
1416           s += sprintf("%4d %-7s %-4s %-4s %-2s %7.3f %7.3f %7.3f %6.3f [%s]\n",
1417                 ap.index, sprintf("%3s.%d", ap.res_name, ap.res_seq),
1418                 quote(ap.name), quote(ap.atom_type), quote(ap.element),
1419                 ap.r.x, ap.r.y, ap.r.z, ap.charge,
1420                 ap.connects.join(","))
1421         }
1422         print s
1423   end
1424
1425   def from_dump(arg)
1426     if natoms > 0
1427           raise "Molecule must be empty"
1428         end
1429     format = "index residue name atom_type element rx ry rz charge connects"
1430         keys = []
1431         resAtoms = Hash.new
1432         newBonds = []
1433         #  arg can be either a String or an array of String.
1434         #  Iterates for each line in the string or each member of the array.
1435         if arg.is_a?(String)
1436           arg = arg.split("\n")
1437         end
1438     arg.each { |line|
1439           if line =~ /^\#/
1440             format = line[1..-1]
1441                 keys = []
1442           end
1443           if keys.length == 0
1444             keys = format.split(" ").collect { |key| key.to_sym }
1445           end
1446           values = line.chomp.split(" ")
1447           next if values == nil || values.length == 0
1448           ap = create_atom(sprintf("X%03d", natoms))
1449           r = Vector3D[0, 0, 0]
1450           keys.each_index { |i|
1451             break if (value = values[i]) == nil
1452                 if value == "\"\""
1453                   value = ""
1454                 else
1455                   value.gsub(/%09/, "\t")
1456                   value.gsub(/%20/, " ")
1457                   value.gsub(/%%/, "%")
1458                 end
1459                 key = keys[i]
1460                 if key == :residue
1461                   if resAtoms[value] == nil
1462                     resAtoms[value] = []
1463                   end
1464                   resAtoms[value].push(ap.index)
1465                 elsif key == :rx
1466                   r.x = value.to_f
1467                 elsif key == :ry
1468                   r.y = value.to_f
1469                 elsif key == :rz
1470                   r.z = value.to_f
1471                 elsif key == :connects
1472                   value.scan(/\d+/).each { |i|
1473                     i = i.to_i
1474                     if ap.index < i
1475                           newBonds.push(ap.index)
1476                           newBonds.push(i)
1477                         end
1478                   }
1479                 elsif key == :index
1480                   next
1481                 else
1482                   ap.set_attr(key, value)
1483                 end
1484           }
1485           ap.r = r
1486         }
1487         resAtoms.each_key { |key|
1488           assign_residue(atom_group(resAtoms[key]), key)
1489         }
1490         create_bond(*newBonds)
1491         self
1492   end
1493
1494   #  Plug-in for loading mbsf
1495   def loadmbsf_plugin(s, lineno)
1496     ""
1497   end
1498
1499   #  Plug-in for saving mbsf
1500   def savembsf_plugin
1501     ""
1502   end
1503
1504 end