OSDN Git Service

Molden import: handling of spherical functions is improved
[molby/Molby.git] / Scripts / loadsave.rb
index db3eb00..bdf545c 100755 (executable)
@@ -1,3 +1,4 @@
+# coding: utf-8
 #
 #  loadsave.rb
 #
@@ -143,7 +144,7 @@ class Molecule
                break if line =~ /TOTAL NUMBER OF BASIS SET/
                if line =~ /^\s*$/
                  #  End of one shell
-                 add_gaussian_orbital_shell(sym, nprims, i)
+                 add_gaussian_orbital_shell(i, sym, nprims)
                  nprims = 0
                  sym = -10
                  next
@@ -220,7 +221,7 @@ class Molecule
                        }
                end
                mo.each_with_index { |m, i|
-                       idx = Integer(mo_labels[i]) - 1
+                       idx = Integer(mo_labels[i])
                        set_mo_coefficients(idx + (alpha ? 0 : ncomps), Float(mo_energies[i]), m)
                #       if mo_labels[i] % 8 == 1
                #               puts "set_mo_coefficients #{idx}, #{mo_energies[i]}, [#{m[0]}, ..., #{m[-1]}]"
@@ -264,23 +265,26 @@ class Molecule
                                break
                        end
                        line.chomp!
-                       if line =~ /ATOM\s+ATOMIC\s+COORDINATES/ || line =~ /COORDINATES OF ALL ATOMS ARE/
+                       if line =~ /ATOM\s+ATOMIC\s+COORDINATES/ || line =~ /COORDINATES OF ALL ATOMS ARE/ || line =~ /COORDINATES \(IN ANGSTROM\) FOR \$DATA GROUP ARE/
                                set_progress_message(mes + "\nReading atomic coordinates...")
                                if line =~ /ATOMIC/
                                  first_line = true
-                                 if !new_unit
-                                   next   #  Skip initial atomic coordinates unless loading into an empty molecule
-                                 end
+                               #  if !new_unit
+                               #    next   #  Skip initial atomic coordinates unless loading into an empty molecule
+                               #  end
+                                 line = fp.gets  #  Skip one line
                                else
                                  first_line = false
                                  nsearch += 1
+                                 if line =~ /COORDINATES OF ALL ATOMS ARE/
+                                   line = fp.gets  #  Skip one line
+                             end
                                end
-                               line = fp.gets    #  Skip one line
                                n = 0
                                coords = []
                                names = []
                                while (line = fp.gets) != nil
-                                       break if line =~ /^\s*$/ || line =~ /END OF ONE/
+                                       break if line =~ /^\s*$/ || line =~ /END OF ONE/ || line =~ /\$END/
                                        next if line =~ /-----/
                                        name, charge, x, y, z = line.split
                                        v = Vector3D[x, y, z]
@@ -308,8 +312,21 @@ class Molecule
                                        new_unit = false
                                #       create_frame
                                else
-                                       if search_mode != 1 || nsearch > 1
-                                               #  The first frame for geometry search has the same coordinates as input
+                                   dont_create = false
+                                       if (search_mode == 1 && nsearch == 1) || first_line
+                                               #  The input coordinate and the first frame for geometry search
+                                               #  can have the same coordinate as the last frame; if this is the case, then
+                                               #  do not create the new frame
+                                               select_frame(nframes - 1)
+                                               dont_create = true
+                                               each_atom { |ap|
+                                                 if (ap.r - coords[ap.index]).length2 > 1e-8
+                                                   dont_create = false
+                                                       break
+                                                 end
+                                               }
+                                       end
+                                       if !dont_create
                                                create_frame([coords])  #  Should not be (coords)
                                        end
                                end
@@ -391,13 +408,20 @@ class Molecule
                                end
                                sub_load_gamess_log_mo_coefficients(lines, lineno, ncomps)
                                set_progress_message(mes)
+                       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/
+                               nbo_lines = []
+                               while (line = fp.gets) != nil
+                                 break if line =~ /done with NBO analysis/
+                                 nbo_lines.push(line)
+                               end
+                               import_nbo_log(nbo_lines)
                        end
                end
        end
        if nframes > 0
          select_frame(nframes - 1)
        end
-       if energy != 0.0
+       if energy && energy != 0.0
          set_property("energy", energy)
        end
        hide_progress_panel
@@ -487,18 +511,590 @@ class Molecule
        (n > 0 ? true : false)
   end
 
+  def sub_load_psi4_log(fp)
+    if natoms == 0
+      new_unit = true
+    else
+      new_unit = false
+    end
+    n = 0
+    nf = 0
+    energy = nil
+  
+    show_progress_panel("Loading Psi4 output file...")
+
+    getline = lambda { @lineno += 1; return fp.gets }
+
+    #  Import coordinates and energies
+    vecs = []
+    ats = []
+    first_frame = nframes
+    trans = nil
+    hf_type = nil
+    nalpha = nil
+    nbeta = nil
+    while line = getline.call
+      if line =~ /==> Geometry <==/
+        #  Skip until line containing "------"
+        while line = getline.call
+          break if line =~ /------/
+        end
+        vecs.clear
+        index = 0
+        #  Read atom positions
+        while line = getline.call
+          line.chomp!
+          break if line =~ /^\s*$/
+          tokens = line.split(' ')
+          if natoms > 0 && first_frame == nframes
+            if index >= natoms || tokens[0] != atoms[index].element
+              hide_progress_panel
+              raise MolbyError, "The atom list does not match the current structure at line #{@lineno}"
+            end
+          end
+          vecs.push(Vector3D[Float(tokens[1]), Float(tokens[2]), Float(tokens[3])])
+          if natoms == 0
+            ats.push(tokens[0])
+          end
+          index += 1
+        end
+        if natoms == 0
+          #  Create molecule from the initial geometry
+          ats.each_with_index { |aname, i|
+            #  Create atoms
+            ap = add_atom(aname)
+            ap.element = aname
+            ap.atom_type = ap.element
+            ap.name = sprintf("%s%d", aname, i)
+            ap.r = vecs[i]
+          }
+          guess_bonds
+        else
+          if vecs.length != natoms
+            break  #  Log file is incomplete
+          end
+          #  Does this geometry differ from the last one?
+          vecs.length.times { |i|
+            if (atoms[i].r - vecs[i]).length2 > 1.0e-14
+              #  Create a new frame and break
+              create_frame
+              vecs.length.times { |j|
+                atoms[j].r = vecs[j]
+              }
+              break
+            end
+          }
+        end
+        #  end geometry
+      elsif line =~ /Final Energy: +([-.0-9]+)/
+        #  Energy for this geometry
+        energy = Float($1)
+        set_property("energy", energy)
+        if line =~ /RHF/
+          hf_type = "RHF"
+        elsif line =~ /UHF/
+          hf_type = "UHF"
+        elsif line =~ /ROHF/
+          hf_type = "ROHF"
+        end
+      elsif line =~ /^ *Nalpha *= *(\d+)/
+        nalpha = Integer($1)
+      elsif line =~ /^ *Nbeta *= *(\d+)/
+        nbeta = Integer($1)
+      end
+    end
+    hide_progress_panel
+    clear_basis_set
+    clear_mo_coefficients
+    set_mo_info(:type => hf_type, :alpha => nalpha, :beta => nbeta)
+    return true
+  end
+
+  #  mol.set_mo_info should be set before calling this function
+  #  Optional label is for importing JANPA output: "NAO" or "CPLO"
+  #  If label is not nil, then returns a hash containing the following key/value pairs:
+  #    :atoms => an array of [element_symbol, seq_num, atomic_num, x, y, z] (angstrom)
+  #    :gto => an array of an array of [sym, [ex0, c0, ex1, c1, ...]]
+  #    :moinfo => an array of [sym, energy, spin (0 or 1), occ]
+  #    :mo => an array of [c0, c1, ...]
+  def sub_load_molden(fp, label = nil)
+    getline = lambda { @lineno += 1; return fp.gets }
+    bohr = 0.529177210903
+    errmsg = nil
+    ncomps = 0  #  Number of components (AOs)
+    occ_alpha = 0  #  Number of occupied alpha orbitals
+    occ_beta = 0   #  Number of occupied beta orbitals
+    if label
+      hash = Hash.new
+    end
+    #  The GTOs (orbital type, contractions and exponents) are stored in gtos[]
+    #  and set just before first [MO] is processed.
+    #  This is because we do not know whether the orbital type is cartesian or spherical
+    #  until we see lines like "[5D]".
+    gtos = []
+    spherical_d = false
+    spherical_f = false
+    spherical_g = false
+    #  Number of components for each orbital type
+    ncomp_hash = { 0=>1, 1=>3, -1=>4, 2=>6, -2=>5, 3=>10, -3=>7, 4=>15, -4=>9 }
+    catch :ignore do
+      while line = getline.call
+        if line =~ /^\[Atoms\]/
+          i = 0
+          while line = getline.call
+            if line =~ /^[A-Z]/
+              #  element, index, atomic_number, x, y, z (in AU)
+              a = line.split(' ')
+              if label
+                (hash[:atoms] ||= []).push([a[0], Integer(a[1]), Integer(a[2]), Float(a[3]) * bohr, Float(a[4]) * bohr, Float(a[5]) * bohr])
+              else
+                if atoms[i].atomic_number != Integer(a[2]) ||
+                  (atoms[i].x - Float(a[3]) * bohr).abs > 1e-4 ||
+                  (atoms[i].y - Float(a[4]) * bohr).abs > 1e-4 ||
+                  (atoms[i].z - Float(a[5]) * bohr).abs > 1e-4
+                  errmsg = "The atom list does not match the current molecule."
+                  throw :ignore
+                end
+              end
+              i += 1
+            else
+              break
+            end
+          end
+          redo  #  The next line will be the beginning of the next block
+        elsif line =~ /^\[GTO\]/
+          shell = 0
+          atom_index = 0
+          while line = getline.call
+            #  index, 0?
+            a = line.split(' ')
+            break if a.length != 2
+            atom_gtos = []  #  [[sym1, [e11, c11, e12, c12, ...], add_exp1], [sym2, [e21, c22, ...], add_exp2], ...]
+            #  loop for shells
+            while line = getline.call
+              #  type, no_of_primitives, 1.00?
+              a = line.split(' ')
+              break if a.length != 3   #  Terminated by a blank line
+              a[0] =~ /^([a-z]+)([0-9]+)?$/
+              symcode = $1
+              add_exp = ($2 == nil ? 0 : $2.to_i)
+              case symcode
+              when "s"
+                sym = 0
+              when "p"
+                sym = 1
+              when "d"
+                sym = 2
+              when "f"
+                sym = 3
+              when "g"
+                sym = 4
+              else
+                raise MolbyError, "Unknown gaussian shell type '#{a[0]}' at line #{@lineno} in MOLDEN file"
+              end
+              nprimitives = Integer(a[1])
+              gtoline = [sym, [], add_exp]
+              atom_gtos.push(gtoline)
+              nprimitives.times { |i|
+                line = getline.call   #  exponent, contraction
+                b = line.split(' ')
+                gtoline[1].push(Float(b[0]), Float(b[1]))
+              }
+              #  end of one shell
+              shell += 1
+            end
+            #  end of one atom
+            atom_index += 1
+            gtos.push(atom_gtos)
+          end
+          if label
+            hash[:gto] = gtos
+          end
+          redo  #  The next line will be the beginning of the next block
+        elsif line =~ /^\[5D\]/ || line =~ /^\[5D7F\]/
+          spherical_d = spherical_f = true
+        elsif line =~ /^\[5D10F\]/
+          spherical_d = true
+          spherical_f = false
+        elsif line =~ /^\[7F\]/
+          spherical_f = true
+        elsif line =~ /^\[9G\]/
+          spherical_g = true
+        elsif line =~ /^\[MO\]/
+          #  Add shell info and primitive coefficients to molecule
+          gtos.each_with_index { | atom_gtos, atom_index|
+            atom_gtos.each { |gtoline|
+              sym = gtoline[0]
+              #  Change orbital type if we use spherical functions
+              sym = -2 if sym == 2 && spherical_d
+              sym = -3 if sym == 3 && spherical_f
+              sym = -4 if sym == 4 && spherical_g
+              gtoline[0] = sym
+              coeffs = gtoline[1]
+              nprimitives = coeffs.length / 2
+              add_exp = gtoline[2]
+              ncomps += ncomp_hash[sym]
+              if !label
+                add_gaussian_orbital_shell(atom_index, sym, nprimitives, add_exp)
+                nprimitives.times { |prim|
+                  add_gaussian_primitive_coefficients(coeffs[prim * 2], coeffs[prim * 2 + 1], 0.0)
+                }
+              end
+            }
+          }
+          m = []
+          idx_alpha = 1   #  set_mo_coefficients() accepts 1-based index of MO
+          idx_beta = 1
+          if label
+            hash[:mo] = []
+            hash[:moinfo] = []
+          end
+          while true
+            #  Loop for each MO
+            m.clear
+            ene = nil
+            spin = nil
+            sym = nil   #  Not used in Molby
+            occ = nil
+            i = 0
+            while line = getline.call
+              if line =~ /^ *Sym= *(\w+)/
+                sym = $1
+              elsif line =~ /^ *Ene= *([-+.0-9e]+)/
+                ene = Float($1)
+              elsif line =~ /^ *Spin= *(\w+)/
+                spin = $1
+              elsif line =~ /^ *Occup= *([-+.0-9e]+)/
+                occ = Float($1)
+                if occ > 0.0
+                  if spin == "Alpha"
+                    occ_alpha += 1
+                  else
+                    occ_beta += 1
+                  end
+                end
+                if label
+                  hash[:moinfo].push([sym, ene, (spin == "Alpha" ? 0 : 1), occ])
+                end
+              elsif line =~ /^ *([0-9]+) +([-+.0-9e]+)/
+                m[i] = Float($2)
+                i += 1
+                if i >= ncomps
+                  if spin == "Alpha"
+                    idx = idx_alpha
+                    idx_alpha += 1
+                  else
+                    idx = idx_beta
+                    idx_beta += 1
+                  end
+                  set_mo_coefficients(idx, ene, m)
+                  if label
+                    hash[:mo].push(m.dup)
+                  end
+                  break
+                end
+              else
+                break
+              end
+            end
+            break if i < ncomps  #  no MO info was found
+          end
+          #  TODO: reorder D, F, G coefficients for Molby order
+          next
+        end #  end if
+      end   #  end while
+    end     #  end catch
+    if errmsg
+      message_box("The MOLDEN file was found but not imported. " + errmsg, "Psi4 import info", :ok)
+      return (label ? nil : false)
+    end
+    return (label ? hash : true)
+  end
+
+  #  Import the JANPA log and related molden files
+  #  Files: inppath.{NAO.molden,CLPO.molden,janpa.log}
+  #  If inppath.spherical.molden is available, then clear existing mo info
+  #  and load from it (i.e. use the basis set converted by molden2molden)
+  def sub_load_janpa_log(inppath)
+    begin
+      fp = File.open(inppath + ".janpa.log", "rt") rescue fp = nil
+      if fp == nil
+        hide_progress_panel  #  Close if it is open
+        message_box("Cannot open JANPA log file #{inppath + '.janpa.log'}: " + $!.to_s)
+        return false
+      end
+      print("Importing #{inppath}.janpa.log.\n")
+      lineno = 0
+      getline = lambda { lineno += 1; return fp.gets }
+      h = Hash.new
+      mfiles = Hash.new
+      h["software"] = "JANPA"
+      nao_num = nil  #  Set later
+      nao_infos = [] #  index=atom_index, value=Hash with key "s", "px", "py" etc.
+      #  nao_infos[index][key]: array of [nao_num, occupancy], in the reverse order of appearance
+      while line = getline.call
+        if line =~ /molden2molden: a conversion tool for MOLDEN/
+          while line = getline.call
+            break if line =~ /^All done!/
+            if line =~ /\.spherical\.molden/
+              #  The MOs are converted to spherical basis set
+              #  Clear the existing MO and load *.spherical.molden
+              sname = inppath + ".spherical.molden"
+              fps = File.open(sname, "rt") rescue fps = nil
+              if fps != nil
+                print("Importing #{sname}.\n")
+                @lineno = 0
+                type = get_mo_info(:type)
+                alpha = get_mo_info(:alpha)
+                beta = get_mo_info(:beta)
+                clear_basis_set
+                set_mo_info(:type=>type, :alpha=>alpha, :beta=>beta)
+                #  mol.@hf_type should be set before calling sub_load_molden
+                @hf_type = type
+                sub_load_molden(fps)
+                fps.close
+              end
+            end
+          end
+        elsif line =~ /^NAO \#/
+          h["NAO"] = []
+          while line = getline.call
+            break if line !~ /^\s*[1-9]/
+            num = Integer(line[0, 5])
+            name = line[5, 21]
+            occ = Float(line[26, 11])
+            #  like A1*: R1*s(0)
+            #  atom_number, occupied?, group_number, orb_sym, angular_number
+            name =~ /\s*[A-Z]+([0-9]+)(\*?):\s* R([0-9]+)\*([a-z]+)\(([-0-9]+)\)/
+            anum = Integer($1)
+            occupied = $2
+            group_num = Integer($3)
+            orb_sym = $4
+            ang_num = Integer($5)
+            orb_desc = orb_sym
+            if orb_desc == "p"
+              orb_desc += ["z", "x", "y"][ang_num + 1]
+            elsif orb_desc == "d"
+            #  TODO: handle d, f, g orbitals
+            end
+            h["NAO"].push([num, anum, occupied, group_num, orb_desc, occ])
+            nao_num = h["NAO"].length
+            ((nao_infos[anum - 1] ||= Hash.new)[orb_desc] ||= []).unshift([nao_num, occ])
+          end
+          #  Create labels
+          h["NAO_L"] = []
+          nao_infos.each_with_index { |value, atom_index|
+            aname = self.atoms[atom_index].name
+            value.each { |orb_desc, ar|
+              ar.each_with_index { |v, group_index|
+                if v[1] > 1.9
+                  label = "core"
+                elsif v[1] > 0.01
+                  label = "val"
+                else
+                  label = "ryd"
+                end
+                principle = group_index + 1
+                orb_sym = orb_desc[0]
+                if orb_sym == "p"
+                  principle += 1
+                elsif orb_sym == "d"
+                  principle += 2
+                elsif orb_sym == "f"
+                  principle += 3
+                elsif orb_sym == "g"
+                  principle += 4
+                end
+                h["NAO_L"][v[0] - 1] = "#{aname} (#{principle}#{orb_desc}) (#{label})"
+              }
+            }
+          }
+        elsif line =~ /^\s*(C?)LPO\s+D e s c r i p t i o n\s+Occupancy\s+Composition/
+          if $1 == "C"
+            key = "CLPO"
+          else
+            key = "LPO"
+          end
+          h[key] = []
+          while line = getline.call
+            break if line =~ /^\s*$/
+            num = Integer(line[0, 5])
+            label1 = line[5, 6].strip
+            desc = line[11, 30].strip
+            occ = line[41, 11].strip
+            comp = line[52, 1000].strip
+            desc =~ /\s*([-A-Za-z0-9]+)(,\s*(.*$))?/
+            desc1 = $1
+            desc2 = ($3 || "")
+            if desc2 =~ /^(.*)*\(NB\)\s*$/ && label1 == ""
+              label1 = "(NB)"
+              desc2 = $1.strip
+            end
+            atoms = desc1.scan(/[A-Za-z]+(\d+)/)   # "C1-H3" -> [["1"], ["3"]]
+            atoms = atoms.map { |a| Integer(a[0]) }  # [1, 3]
+            hybrids_a = comp.scan(/h(\d+)@[A-Za-z]+(\d+)/)  #  "h8@C1...h13@H3" -> "[["8", "1"], ["13", "3"]]
+            hybrids = []
+            hybrids_a.each { |a|
+              i = atoms.find_index(Integer(a[1]))
+              if i != nil
+                hybrids[i] = Integer(a[0])
+              end
+            } # [8, 13]
+            #  like ["(BD)", [1, 3], "Io = 0.2237", occ, [8, 13]]
+            #  1, 3 are the atom indices (1-based)
+            #  8, 13 are the number of hybrid orbitals (1-based)
+            h[key][num - 1] = [label1, atoms, desc2, Float(occ), hybrids]
+          end
+          h[key + "_L"] = []
+          if key == "CLPO"
+            #  Also register labels of "LHO"
+            h["LHO_L"] = [""] * nao_num
+          end
+          nao_num.times { |i|
+            val = h[key][i]
+            if val == nil
+              label = ""  #  The labels for Rydberg orbitals may be replaced later
+            else
+              aname1 = self.atoms[val[1][0] - 1].name rescue aname1 = ""
+              aname2 = self.atoms[val[1][1] - 1].name rescue aname2 = ""
+              if aname2 == ""
+                label = "#{aname1} #{val[0]}"
+              else
+                label = "#{aname1}(#{aname2}) #{val[0]}"
+              end
+            end
+            h[key + "_L"][i] = label
+            if key == "CLPO" && val != nil && val[0] != "(NB)"
+              hybrids = val[4]
+              kind = (val[0] == "(BD)" ? "(val)" : "(lp)")
+              if aname2 == ""
+                label = "#{aname1} #{kind}"
+              else
+                label = "#{aname1}(#{aname2}) #{kind}"
+              end
+              h["LHO_L"][hybrids[0] - 1] = label
+              if hybrids[1] != nil
+                #  aname2 should be non-empty
+                label = "#{aname2}(#{aname1}) #{kind}"
+                h["LHO_L"][hybrids[1] - 1] = label
+              end
+            end
+          }
+        elsif line =~ /^ -NAO_Molden_File: (\S*)/
+          mfiles["NAO"] = $1
+        elsif line =~ /^ -LHO_Molden_File: (\S*)/
+          mfiles["LHO"] = $1
+        elsif line =~ /^ -CLPO_Molden_File: (\S*)/
+          mfiles["CLPO"] = $1
+        elsif line =~ /^ -PNAO_Molden_File: (\S*)/
+          mfiles["PNAO"] = $1
+        elsif line =~ /^ -AHO_Molden_File: (\S*)/
+          mfiles["AHO"] = $1
+        elsif line =~ /^ -LPO_Molden_File: (\S*)/
+          mfiles["LPO"] = $1
+        end
+      end
+      fp.close
+      #  Read molden files
+      mfiles.each { |key, value|
+        fp = Kernel.open(value, "rt") rescue fp = nil
+        if fp
+          print("Importing #{value}.\n")
+          res = sub_load_molden(fp, key)
+          if res
+            #  Some kind of orbital based on AO
+            h["AO/#{key}"] = LAMatrix.new(res[:mo])
+          end
+          fp.close
+          if key == "CLPO" || key == "LPO" || key == "LHO"
+            #  Set the label of Rydberg orbitals if possible
+            if h[key + "_L"] != nil
+              a = h["AO/#{key}"]
+              nao_num.times { |i|
+                label = h[key + "_L"][i]
+                if label == ""
+                  max_idx = nil
+                  max_val = -1.0
+                  nao_infos.each_with_index { |inf, atom_index|
+                    atomic_contrib = 0.0
+                    inf.each { |k, v| # k is "s", "px" etc, v is array of [nao_num, occupancy]
+                      #  Sum for all naos belonging to this atom
+                      v.each { |num_occ|
+                        atomic_contrib += a[i, num_occ[0] - 1] ** 2
+                      }
+                    }
+                    if atomic_contrib > max_val
+                      max_val = atomic_contrib
+                      max_idx = atom_index
+                    end
+                  }
+                  label = self.atoms[max_idx].name + " (ry)"
+                  h[key + "_L"][i] = label
+                end
+              }
+            end
+          end
+        end
+      }
+      @nbo = h
+      if @nbo["AO/NAO"] && @nbo["AO/LHO"] && @nbo["AO/PNAO"]
+        #  Generate PLHO from PNAO, NAO, LHO
+        #  This protocol was suggested by the JANPA author in a private commnunication.
+        begin
+          nao2lho = @nbo["AO/NAO"].inverse * @nbo["AO/LHO"]
+          nao2pnao = @nbo["AO/NAO"].inverse * @nbo["AO/PNAO"]
+          sign = LAMatrix.diagonal((0...nao2pnao.column_size).map { |i| (nao2pnao[i, i] < 0 ? -1 : 1)})
+          @nbo["AO/PLHO"] = @nbo["AO/PNAO"] * sign * nao2lho
+        rescue
+          @nbo["AO/PLHO"] = nil
+        end
+      end
+      return true
+    rescue => e
+      $stderr.write(e.message + "\n")
+      $stderr.write(e.backtrace.inspect + "\n")
+    end
+  end
+
   def loadout(filename)
-    retval = false
-    fp = open(filename, "rb")
-       while s = fp.gets
-         if s =~ /Gaussian/
-           retval = sub_load_gaussian_log(fp)
-               break
-         elsif s =~ /GAMESS/
-           retval = sub_load_gamess_log(fp)
-               break
-         end
-       end
+  retval = false
+  fp = open(filename, "rb")
+  @lineno = 0
+  begin
+    while s = fp.gets
+      @lineno += 1
+      if s =~ /Gaussian/
+        retval = sub_load_gaussian_log(fp)
+        break
+      elsif s =~ /GAMESS/
+        retval = sub_load_gamess_log(fp)
+        break
+      elsif s =~ /Psi4/
+        retval = sub_load_psi4_log(fp)
+        if retval
+          #  If .molden file exists, then try to read it
+          namepath = filename.gsub(/\.\w*$/, "")
+          mname = "#{namepath}.molden"
+          if File.exists?(mname)
+            fp2 = open(mname, "rb")
+            if fp2
+              flag = sub_load_molden(fp2)
+              fp2.close
+              status = (flag ? 0 : -1)
+            end
+          end
+          if File.exists?("#{namepath}.janpa.log")
+            flag = sub_load_janpa_log(namepath)
+            status = (flag ? 0 : -1)
+          end
+        end
+        break
+      end
+    end
+  rescue
+    hide_progress_panel
+    raise
+  end
        fp.close
        return retval
   end
@@ -738,6 +1334,7 @@ end_of_header
   alias :savegjf :savecom
   
   def loadcif(filename)
+    mol = self
     def getciftoken(fp)
          while @tokens.length == 0
            line = fp.gets
@@ -797,302 +1394,329 @@ end_of_header
          end
          sym
        end
-       def find_atom_by_name(name)
+       def find_atom_by_name(mol, name)
          name = name.delete(" ()")
-         ap = self.atoms[name] rescue ap = nil
+         ap = mol.atoms[name] rescue ap = nil
          return ap
        end
-       warn_message = ""
-       verbose = nil
-       bond_defined = false
-       @tokens = []
-       special_positions = []
-       self.remove(All)
+       selfname = self.name
        fp = open(filename, "rb")
-       cell = []
-       cell_trans = cell_trans_inv = Transform.identity
-       token = getciftoken(fp)
-       pardigits_re = /\(\d+\)/
-       calculated_atoms = []
-       while token != nil
-         if token =~ /^_cell/
-               val = getciftoken(fp)
-               if token == "_cell_length_a"
-                 cell[0], cell[6] = float_strip_rms(val)
-               elsif token == "_cell_length_b"
-                 cell[1], cell[7] = float_strip_rms(val)
-               elsif token == "_cell_length_c"
-                 cell[2], cell[8] = float_strip_rms(val)
-               elsif token == "_cell_angle_alpha"
-                 cell[3], cell[9] = float_strip_rms(val)
-               elsif token == "_cell_angle_beta"
-                 cell[4], cell[10] = float_strip_rms(val)
-               elsif token == "_cell_angle_gamma"
-                 cell[5], cell[11] = float_strip_rms(val)
-               end
-               if cell.length == 12 && cell.all?
-                 self.cell = cell
-                 puts "Unit cell is set to #{cell.inspect}." if verbose
-                 cell = []
-                 cell_trans = self.cell_transform
-                 cell_trans_inv = cell_trans.inverse
-               end
-               token = getciftoken(fp)
-               next
-      elsif token.casecmp("#loop_") == 0
-           labels = []
-               while (token = getciftoken(fp)) && token[0] == ?_
-                 labels.push(token)
-               end
-               if labels[0] =~ /symmetry_equiv_pos|atom_site_label|atom_site_aniso_label|geom_bond/
-                 hlabel = Hash.new(-10000000)
-                 labels.each_with_index { |lb, i|
-                       hlabel[lb] = i
-                 }
-                 data = []
-                 n = labels.length
-                 a = []
-                 while 1
-                       break if token == nil || token[0] == ?_ || token[0] == ?#
-                       a.push(token)
-                       if a.length == n
-                         data.push(a)
-                         a = []
-                       end
+       data_identifier = nil
+       @tokens = []
+       count_up = 1
+       while true
+         warn_message = ""
+         verbose = nil
+         bond_defined = false
+         special_positions = []
+         mol.remove(All)
+         cell = []
+         cell_trans = cell_trans_inv = Transform.identity
+         token = getciftoken(fp)
+         pardigits_re = /\(\d+\)/
+         calculated_atoms = []
+         while token != nil
+           if token =~ /^\#data_/i
+                 if data_identifier == nil || mol.natoms == 0
+                   #  First block or no atoms yet
+            #  Continue processing of this molecule
+                   data_identifier = token
                        token = getciftoken(fp)
+                       next
+                 else
+                   #  Description of another molecule begins here
+                       data_identifier = token
+                       break
+                 end
+           elsif token =~ /^_cell/
+                 val = getciftoken(fp)
+                 if token == "_cell_length_a"
+                   cell[0], cell[6] = float_strip_rms(val)
+                 elsif token == "_cell_length_b"
+                   cell[1], cell[7] = float_strip_rms(val)
+                 elsif token == "_cell_length_c"
+                   cell[2], cell[8] = float_strip_rms(val)
+                 elsif token == "_cell_angle_alpha"
+                   cell[3], cell[9] = float_strip_rms(val)
+                 elsif token == "_cell_angle_beta"
+                   cell[4], cell[10] = float_strip_rms(val)
+                 elsif token == "_cell_angle_gamma"
+                   cell[5], cell[11] = float_strip_rms(val)
                  end
-                 if labels[0] =~ /^_symmetry_equiv_pos/
-                   data.each { |d|
-                         symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]]
-                         symstr.delete("\"\'")
-                         exps = symstr.split(/,/)
-                         sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
-                         exps.each_with_index { |s, i|
-                           terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
-                               terms.each { |a|
-                                 #  a[0]: sign, a[2]: numerator, a[4]: denometer
-                                 if a[4] != nil
-                                   #  The number part is a[2]/a[4]
-                                   num = Float(a[2])/Float(a[4])
-                                 elsif a[2] != nil
-                                   #  The number part is either integer or a floating point
-                                   num = Float(a[2])
+                 if cell.length == 12 && cell.all?
+                   mol.cell = cell
+            puts "Unit cell is set to #{mol.cell.inspect}." if verbose
+                   cell = []
+                   cell_trans = mol.cell_transform
+                   cell_trans_inv = cell_trans.inverse
+                 end
+                 token = getciftoken(fp)
+                 next
+        elsif token.casecmp("#loop_") == 0
+             labels = []
+                 while (token = getciftoken(fp)) && token[0] == ?_
+                   labels.push(token)
+                 end
+                 if labels[0] =~ /symmetry_equiv_pos|space_group_symop|atom_site_label|atom_site_aniso_label|geom_bond/
+                   hlabel = Hash.new(-10000000)
+                   labels.each_with_index { |lb, i|
+                         hlabel[lb] = i
+                   }
+                   data = []
+                   n = labels.length
+                   a = []
+                   while true
+                         break if token == nil || token[0] == ?_ || token[0] == ?#
+                         a.push(token)
+                         if a.length == n
+                           data.push(a)
+                           a = []
+                         end
+                         token = getciftoken(fp)
+                   end
+                   if labels[0] =~ /^_symmetry_equiv_pos/ || labels[0] =~ /^_space_group_symop/
+                     data.each { |d|
+                           symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]] || d[hlabel["_space_group_symop_operation_xyz"]]
+                           symstr.delete("\"\'")
+                           exps = symstr.split(/,/)
+                           sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
+                           exps.each_with_index { |s, i|
+                             terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
+                                 terms.each { |a|
+                                   #  a[0]: sign, a[2]: numerator, a[4]: denometer
+                                   if a[4] != nil
+                                     #  The number part is a[2]/a[4]
+                                     num = Float(a[2])/Float(a[4])
+                                   elsif a[2] != nil
+                                     #  The number part is either integer or a floating point
+                                     num = Float(a[2])
+                                   else
+                                     num = 1.0
+                                   end
+                                   num = -num if a[0][0] == ?-
+                                   xyz = (a[5] || a[6])
+                                   if xyz == "x" || xyz == "X"
+                                     sym[i] = num
+                                   elsif xyz == "y" || xyz == "Y"
+                                     sym[i + 3] = num
+                                   elsif xyz == "z" || xyz == "Z"
+                                     sym[i + 6] = num
+                                   else
+                                     sym[9 + i] = num
+                                   end
+                                 }
+                           }
+                           puts "symmetry operation #{sym.inspect}" if verbose
+                           mol.add_symmetry(Transform.new(sym))
+                         }
+                         puts "#{mol.nsymmetries} symmetry operations are added" if verbose
+                   elsif labels[0] =~ /^_atom_site_label/
+                         #  Create atoms
+                         data.each { |d|
+                           name = d[hlabel["_atom_site_label"]]
+                           elem = d[hlabel["_atom_site_type_symbol"]]
+                           fx = d[hlabel["_atom_site_fract_x"]]
+                           fy = d[hlabel["_atom_site_fract_y"]]
+                           fz = d[hlabel["_atom_site_fract_z"]]
+                           uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
+                           biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
+                           occ = d[hlabel["_atom_site_occupancy"]]
+                           calc = d[hlabel["_atom_site_calc_flag"]]
+                           name = name.delete(" ()")
+                           if elem == nil || elem == ""
+                             if name =~ /[A-Za-z]{1,2}/
+                                   elem = $&.capitalize
                                  else
-                                   num = 1.0
+                                   elem = "Du"
                                  end
-                                 num = -num if a[0][0] == ?-
-                                 xyz = (a[5] || a[6])
-                                 if xyz == "x" || xyz == "X"
-                                   sym[i] = num
-                                 elsif xyz == "y" || xyz == "Y"
-                                   sym[i + 3] = num
-                                 elsif xyz == "z" || xyz == "Z"
-                                   sym[i + 6] = num
-                                 else
-                                   sym[9 + i] = num
+                           end
+                           ap = mol.add_atom(name, elem, elem)
+                           ap.fract_x, ap.sigma_x = float_strip_rms(fx)
+                           ap.fract_y, ap.sigma_y = float_strip_rms(fy)
+                           ap.fract_z, ap.sigma_z = float_strip_rms(fz)
+                           if biso
+                             ap.temp_factor, sig = float_strip_rms(biso)
+                           elsif uiso
+                             ap.temp_factor, sig = float_strip_rms(uiso)
+                                 ap.temp_factor *= 78.9568352087149          #  8*pi*pi
+                           end
+                           ap.occupancy, sig = float_strip_rms(occ)
+                           if calc == "c" || calc == "calc"
+                             calculated_atoms.push(ap.index)
+                       end
+                           #  Guess special positions
+                           (1...mol.nsymmetries).each { |isym|
+                             sr = ap.fract_r
+                             sr = (mol.transform_for_symop(isym) * sr) - sr;
+                                 nx = (sr.x + 0.5).floor
+                                 ny = (sr.y + 0.5).floor
+                                 nz = (sr.z + 0.5).floor
+                                 if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
+                                   #  [isym, -nx, -ny, -nz] transforms this atom to itself
+                                   #  The following line is equivalent to:
+                                   #    if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
+                                   #    special_positions[ap.index].push(...)
+                                   (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
                                  end
-                               }
+                           }
+                           if verbose && special_positions[ap.index]
+                             puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
+                           end
                          }
-                         puts "symmetry operation #{sym.inspect}" if verbose
-                         add_symmetry(Transform.new(sym))
-                       }
-                       puts "#{self.nsymmetries} symmetry operations are added" if verbose
-                 elsif labels[0] =~ /^_atom_site_label/
-                       #  Create atoms
-                       data.each { |d|
-                         name = d[hlabel["_atom_site_label"]]
-                         elem = d[hlabel["_atom_site_type_symbol"]]
-                         fx = d[hlabel["_atom_site_fract_x"]]
-                         fy = d[hlabel["_atom_site_fract_y"]]
-                         fz = d[hlabel["_atom_site_fract_z"]]
-                         uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
-                         biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
-                         occ = d[hlabel["_atom_site_occupancy"]]
-                         calc = d[hlabel["_atom_site_calc_flag"]]
-                         name = name.delete(" ()")
-                         if elem == nil || elem == ""
-                           if name =~ /[A-Za-z]{1,2}/
-                                 elem = $&.capitalize
-                               else
-                                 elem = "Du"
-                               end
-                         end
-                         ap = self.add_atom(name, elem, elem)
-                         ap.fract_x, ap.sigma_x = float_strip_rms(fx)
-                         ap.fract_y, ap.sigma_y = float_strip_rms(fy)
-                         ap.fract_z, ap.sigma_z = float_strip_rms(fz)
-                         if biso
-                           ap.temp_factor, sig = float_strip_rms(biso)
-                         elsif uiso
-                           ap.temp_factor, sig = float_strip_rms(uiso)
-                               ap.temp_factor *= 78.9568352087149          #  8*pi*pi
-                         end
-                         ap.occupancy, sig = float_strip_rms(occ)
-                         if calc == "c" || calc == "calc"
-                           calculated_atoms.push(ap.index)
-                     end
-                         #  Guess special positions
-                         (1...nsymmetries).each { |isym|
-                           sr = ap.fract_r
-                           sr = (transform_for_symop(isym) * sr) - sr;
-                               nx = (sr.x + 0.5).floor
-                               ny = (sr.y + 0.5).floor
-                               nz = (sr.z + 0.5).floor
-                               if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
-                                 #  [isym, -nx, -ny, -nz] transforms this atom to itself
-                                 #  The following line is equivalent to:
-                                 #    if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
-                                 #    special_positions[ap.index].push(...)
-                                 (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
-                               end
+                         puts "#{mol.natoms} atoms are created." if verbose
+                   elsif labels[0] =~ /^_atom_site_aniso_label/
+                     #  Set anisotropic parameters
+                         c = 0
+                         data.each { |d|
+                           name = d[hlabel["_atom_site_aniso_label"]]
+                           ap = find_atom_by_name(mol, name)
+                           next if !ap
+                           u11 = d[hlabel["_atom_site_aniso_U_11"]]
+                           if u11
+                             usig = []
+                             u11, usig[0] = float_strip_rms(u11)
+                             u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
+                             u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
+                             u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
+                             u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
+                             u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
+                             ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
+                                 c += 1
+                           end
                          }
-                         if verbose && special_positions[ap.index]
-                           puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
-                         end
-                       }
-                       puts "#{self.natoms} atoms are created." if verbose
-                 elsif labels[0] =~ /^_atom_site_aniso_label/
-                   #  Set anisotropic parameters
-                       c = 0
-                       data.each { |d|
-                         name = d[hlabel["_atom_site_aniso_label"]]
-                         ap = find_atom_by_name(name)
-                         next if !ap
-                         u11 = d[hlabel["_atom_site_aniso_U_11"]]
-                         if u11
-                           usig = []
-                           u11, usig[0] = float_strip_rms(u11)
-                           u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
-                           u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
-                           u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
-                           u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
-                           u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
-                           ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
-                               c += 1
-                         end
-                       }
-                       puts "#{c} anisotropic parameters are set." if verbose
-                 elsif labels[0] =~ /^_geom_bond/
-                   #  Create bonds
-                       exbonds = []
-                       data.each { |d|
-                         n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
-                         n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
-                         sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
-                         sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
-                         n1 = find_atom_by_name(n1)
-                         n2 = find_atom_by_name(n2)
-                         next if n1 == nil || n2 == nil
-                     n1 = n1.index
-                         n2 = n2.index
-                         sym1 = parse_symmetry_operation(sym1)
-                         sym2 = parse_symmetry_operation(sym2)
-                         if sym1 || sym2
-                           exbonds.push([n1, n2, sym1, sym2])
-                         else
-                           self.create_bond(n1, n2)
-                         end
-                         tr1 = (sym1 ? transform_for_symop(sym1) : Transform.identity)
-                         tr2 = (sym2 ? transform_for_symop(sym2) : Transform.identity)
-                         if special_positions[n1]
-                               #  Add extra bonds for equivalent positions of n1
-                               special_positions[n1].each { |symop|
-                                 sym2x = symop_for_transform(tr1 * transform_for_symop(symop) * tr1.inverse * tr2)
-                                 exbonds.push([n1, n2, sym1, sym2x])
-                               }
-                         end
-                         if special_positions[n2]
-                               #  Add extra bonds n2-n1.symop, where symop transforms n2 to self
-                               tr = (sym1 ? transform_for_symop(sym1) : Transform.identity)
-                               special_positions[n2].each { |symop|
-                                 sym1x = symop_for_transform(tr2 * transform_for_symop(symop) * tr2.inverse * tr1)
-                                 exbonds.push([n2, n1, sym2, sym1x])
-                               }
-                         end                           
-                   }
-                       bond_defined = true
-                       puts "#{self.nbonds} bonds are created." if verbose
-                       if calculated_atoms.length > 0
-                         #  Guess bonds for calculated hydrogen atoms
-                         n1 = 0
-                         calculated_atoms.each { |ai|
-                           if atoms[ai].connects.length == 0
-                                 as = find_close_atoms(ai)
-                                 as.each { |aj|
-                                   self.create_bond(ai, aj)
-                                       n1 += 1
+                         puts "#{c} anisotropic parameters are set." if verbose
+                   elsif labels[0] =~ /^_geom_bond/
+                     #  Create bonds
+                         exbonds = []
+                         data.each { |d|
+                           n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
+                           n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
+                           sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
+                           sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
+                           n1 = find_atom_by_name(mol, n1)
+                           n2 = find_atom_by_name(mol, n2)
+                           next if n1 == nil || n2 == nil
+                       n1 = n1.index
+                           n2 = n2.index
+                           sym1 = parse_symmetry_operation(sym1)
+                           sym2 = parse_symmetry_operation(sym2)
+                           if sym1 || sym2
+                             exbonds.push([n1, n2, sym1, sym2])
+                           else
+                             mol.create_bond(n1, n2)
+                           end
+                           tr1 = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
+                           tr2 = (sym2 ? mol.transform_for_symop(sym2) : Transform.identity)
+                           if special_positions[n1]
+                                 #  Add extra bonds for equivalent positions of n1
+                                 special_positions[n1].each { |symop|
+                                   sym2x = mol.symop_for_transform(tr1 * mol.transform_for_symop(symop) * tr1.inverse * tr2)
+                                   exbonds.push([n1, n2, sym1, sym2x])
                                  }
-                               end
-                         }
-                         puts "#{n1} bonds are guessed." if verbose
-                       end
-                       if exbonds.length > 0
-                         h = Dialog.run("CIF Import: Symmetry Expansion") {
-                           layout(1,
-                                 item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
-                                 item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
-                                 item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
-                                 item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
-                               )
-                               radio_group("atoms_only", "fragment", "ignore")
-                         }
-                         if h[:status] == 0 && h["ignore"] == 0
-                           atoms_only = (h["atoms_only"] != 0)
-                               if !atoms_only
-                                 fragments = []
-                                 self.each_fragment { |f| fragments.push(f) }
-                               end
-                               debug = nil
-                               exbonds.each { |ex|
-                                 if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
-                                 ex0 = ex.dup
-                                 (2..3).each { |i|
-                                   symop = ex[i]
-                                       if symop == nil
-                                         ex[i + 2] = ex[i - 2]
-                                       else
-                                         if debug; puts "  symop = #{symop.inspect}"; end
-                                         #  Expand the atom or the fragment including the atom
-                                         if atoms_only
-                                               ig = IntGroup[ex[i - 2]]
-                                               idx = 0
+                           end
+                           if special_positions[n2]
+                                 #  Add extra bonds n2-n1.symop, where symop transforms n2 to self
+                                 tr = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
+                                 special_positions[n2].each { |symop|
+                                   sym1x = mol.symop_for_transform(tr2 * mol.transform_for_symop(symop) * tr2.inverse * tr1)
+                                   exbonds.push([n2, n1, sym2, sym1x])
+                                 }
+                           end                         
+                     }
+              if mol.nbonds > 0
+                bond_defined = true
+              end
+                         puts "#{mol.nbonds} bonds are created." if verbose
+                         if calculated_atoms.length > 0
+                           #  Guess bonds for calculated hydrogen atoms
+                           n1 = 0
+                           calculated_atoms.each { |ai|
+                             if mol.atoms[ai].connects.length == 0
+                                   as = mol.find_close_atoms(ai)
+                                   as.each { |aj|
+                                     mol.create_bond(ai, aj)
+                                         n1 += 1
+                                   }
+                                 end
+                           }
+                           puts "#{n1} bonds are guessed." if verbose
+                         end
+                         if exbonds.length > 0
+                           h = Dialog.run("CIF Import: Symmetry Expansion") {
+                             layout(1,
+                                   item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
+                                   item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
+                                   item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
+                                   item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
+                                 )
+                                 radio_group("atoms_only", "fragment", "ignore")
+                           }
+                           if h[:status] == 0 && h["ignore"] == 0
+                             atoms_only = (h["atoms_only"] != 0)
+                                 if !atoms_only
+                                   fragments = []
+                                   mol.each_fragment { |f| fragments.push(f) }
+                                 end
+                                 debug = nil
+                                 exbonds.each { |ex|
+                                   if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
+                                   ex0 = ex.dup
+                                   (2..3).each { |i|
+                                     symop = ex[i]
+                                         if symop == nil
+                                           ex[i + 2] = ex[i - 2]
                                          else
-                                               ig = fragments.find { |f| f.include?(ex[i - 2]) }
-                                               ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
+                                           if debug; puts "  symop = #{symop.inspect}"; end
+                                           #  Expand the atom or the fragment including the atom
+                                           if atoms_only
+                                                 ig = IntGroup[ex[i - 2]]
+                                                 idx = 0
+                                           else
+                                                 ig = fragments.find { |f| f.include?(ex[i - 2]) }
+                                                 ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
+                                           end
+                                           symop[4] = ex[i - 2]  #  Base atom
+                                           if debug; puts "  expanding #{ig} by #{symop.inspect}"; end
+                                           a = mol.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
+                                           ex[i + 2] = a[idx]   #  Index of the expanded atom
                                          end
-                                         symop[4] = ex[i - 2]  #  Base atom
-                                         if debug; puts "  expanding #{ig} by #{symop.inspect}"; end
-                                         a = self.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
-                                         ex[i + 2] = a[idx]   #  Index of the expanded atom
-                                       end
+                                   }
+                                   if ex[4] && ex[5] && ex[4] != ex[5]
+                                     if debug; puts "  creating bond #{ex[4]} - #{ex[5]}"; end
+                                     mol.create_bond(ex[4], ex[5])
+                                   end
                                  }
-                                 if ex[4] && ex[5] && ex[4] != ex[5]
-                                   if debug; puts "  creating bond #{ex[4]} - #{ex[5]}"; end
-                                   self.create_bond(ex[4], ex[5])
-                                 end
-                               }
+                           end
                          end
-                       end
-                       puts "#{self.nbonds} bonds are created." if verbose
+                         puts "#{mol.nbonds} bonds are created." if verbose
+                   end
+                   next
+                 else
+                 #  puts "Loop beginning with #{labels[0]} is skipped"
                  end
-                 next
-               else
-               #  puts "Loop beginning with #{labels[0]} is skipped"
-               end
-         else
-           #  Skip this token
-               token = getciftoken(fp)
+           else
+             #  Skip this token
+                 token = getciftoken(fp)
+           end
+           #  Skip tokens until next tag or reserved word is detected
+           while token != nil && token[0] != ?_ && token[0] != ?#
+                 token = getciftoken(fp)
+           end
+           next
+         end
+      if !bond_defined
+               mol.guess_bonds
          end
-         #  Skip tokens until next tag or reserved word is detected
-         while token != nil && token[0] != ?_ && token[0] != ?#
-               token = getciftoken(fp)
+         if token != nil && token == data_identifier
+           #  Process next molecule: open a new molecule and start adding atom on that
+               mol = Molecule.new
+               count_up += 1
+               (@aux_mols ||= []).push(mol)
+               next
          end
-         next
+         break
        end
        fp.close
-       if !bond_defined
-         self.guess_bonds
-       end
 #      self.undo_enabled = save_undo_enabled
        return true
   end
@@ -1129,10 +1753,12 @@ end_of_header
        keys = []
        resAtoms = Hash.new
        newBonds = []
+       #  arg can be either a String or an array of String.
+       #  Iterates for each line in the string or each member of the array.
+       if arg.is_a?(String)
+         arg = arg.split("\n")
+       end
     arg.each { |line|
-      #  arg can be either a String or an array of String. If it is a string,
-         #  arg.each iterates for each line in the string. If it is an array,
-         #  arg.each iterates for each member of the array.
          if line =~ /^\#/
            format = line[1..-1]
                keys = []
@@ -1188,4 +1814,14 @@ end_of_header
        self
   end
 
+  #  Plug-in for loading mbsf
+  def loadmbsf_plugin(s, lineno)
+    ""
+  end
+  
+  #  Plug-in for saving mbsf
+  def savembsf_plugin
+    ""
+  end
+  
 end