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].upcase != atoms[index].element.upcase
 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   #  Optional label is for importing JANPA output: "NAO" or "CPLO"
 615   #  If label is not nil, then returns a hash containing the following key/value pairs:
 616   #    :atoms => an array of [element_symbol, seq_num, atomic_num, x, y, z] (angstrom)
 617   #    :gto => an array of an array of [sym, [ex0, c0, ex1, c1, ...]]
 618   #    :moinfo => an array of [sym, energy, spin (0 or 1), occ]
 619   #    :mo => an array of [c0, c1, ...]
 620   def sub_load_molden(fp, label = nil)
 621     getline = lambda { @lineno += 1; return fp.gets }
 622     bohr = 0.529177210903
 623     errmsg = nil
 624     ncomps = 0  #  Number of components (AOs)
 625     occ_alpha = 0  #  Number of occupied alpha orbitals
 626     occ_beta = 0   #  Number of occupied beta orbitals
 627     if label
 628       hash = Hash.new
 629     end
 630     #  The GTOs (orbital type, contractions and exponents) are stored in gtos[]
 631     #  and set just before first [MO] is processed.
 632     #  This is because we do not know whether the orbital type is cartesian or spherical
 633     #  until we see lines like "[5D]".
 634     gtos = []
 635     spherical_d = false
 636     spherical_f = false
 637     spherical_g = false
 638     #  Number of components for each orbital type
 639     ncomp_hash = { 0=>1, 1=>3, -1=>4, 2=>6, -2=>5, 3=>10, -3=>7, 4=>15, -4=>9 }
 640     catch :ignore do
 641       while line = getline.call
 642         if line =~ /^\[Atoms\]/
 643           i = 0
 644           while line = getline.call
 645             if line =~ /^[A-Z]/
 646               #  element, index, atomic_number, x, y, z (in AU)
 647               a = line.split(' ')
 648               if label
 649                 (hash[:atoms] ||= []).push([a[0], Integer(a[1]), Integer(a[2]), Float(a[3]) * bohr, Float(a[4]) * bohr, Float(a[5]) * bohr])
 650               else
 651                 if atoms[i].atomic_number != Integer(a[2]) ||
 652                   (atoms[i].x - Float(a[3]) * bohr).abs > 1e-4 ||
 653                   (atoms[i].y - Float(a[4]) * bohr).abs > 1e-4 ||
 654                   (atoms[i].z - Float(a[5]) * bohr).abs > 1e-4
 655                   errmsg = "The atom list does not match the current molecule."
 656                   throw :ignore
 657                 end
 658               end
 659               i += 1
 660             else
 661               break
 662             end
 663           end
 664           redo  #  The next line will be the beginning of the next block
 665         elsif line =~ /^\[GTO\]/
 666           shell = 0
 667           atom_index = 0
 668           while line = getline.call
 669             #  index, 0?
 670             a = line.split(' ')
 671             break if a.length != 2
 672             atom_gtos = []  #  [[sym1, [e11, c11, e12, c12, ...], add_exp1], [sym2, [e21, c22, ...], add_exp2], ...]
 673             #  loop for shells
 674             while line = getline.call
 675               #  type, no_of_primitives, 1.00?
 676               a = line.split(' ')
 677               break if a.length != 3   #  Terminated by a blank line
 678               a[0] =~ /^([a-z]+)([0-9]+)?$/
 679               symcode = $1
 680               add_exp = ($2 == nil ? 0 : $2.to_i)
 681               case symcode
 682               when "s"
 683                 sym = 0
 684               when "p"
 685                 sym = 1
 686               when "d"
 687                 sym = 2
 688               when "f"
 689                 sym = 3
 690               when "g"
 691                 sym = 4
 692               else
 693                 raise MolbyError, "Unknown gaussian shell type '#{a[0]}' at line #{@lineno} in MOLDEN file"
 694               end
 695               nprimitives = Integer(a[1])
 696               gtoline = [sym, [], add_exp]
 697               atom_gtos.push(gtoline)
 698               nprimitives.times { |i|
 699                 line = getline.call   #  exponent, contraction
 700                 b = line.split(' ')
 701                 gtoline[1].push(Float(b[0]), Float(b[1]))
 702               }
 703               #  end of one shell
 704               shell += 1
 705             end
 706             #  end of one atom
 707             atom_index += 1
 708             gtos.push(atom_gtos)
 709           end
 710           if label
 711             hash[:gto] = gtos
 712           end
 713           redo  #  The next line will be the beginning of the next block
 714         elsif line =~ /^\[5D\]/ || line =~ /^\[5D7F\]/
 715           spherical_d = spherical_f = true
 716         elsif line =~ /^\[5D10F\]/
 717           spherical_d = true
 718           spherical_f = false
 719         elsif line =~ /^\[7F\]/
 720           spherical_f = true
 721         elsif line =~ /^\[9G\]/
 722           spherical_g = true
 723         elsif line =~ /^\[MO\]/
 724           #  Add shell info and primitive coefficients to molecule
 725           gtos.each_with_index { | atom_gtos, atom_index|
 726             atom_gtos.each { |gtoline|
 727               sym = gtoline[0]
 728               #  Change orbital type if we use spherical functions
 729               sym = -2 if sym == 2 && spherical_d
 730               sym = -3 if sym == 3 && spherical_f
 731               sym = -4 if sym == 4 && spherical_g
 732               gtoline[0] = sym
 733               coeffs = gtoline[1]
 734               nprimitives = coeffs.length / 2
 735               add_exp = gtoline[2]
 736               ncomps += ncomp_hash[sym]
 737               if !label
 738                 add_gaussian_orbital_shell(atom_index, sym, nprimitives, add_exp)
 739                 nprimitives.times { |prim|
 740                   add_gaussian_primitive_coefficients(coeffs[prim * 2], coeffs[prim * 2 + 1], 0.0)
 741                 }
 742               end
 743             }
 744           }
 745           m = []
 746           idx_alpha = 1   #  set_mo_coefficients() accepts 1-based index of MO
 747           idx_beta = 1
 748           if label
 749             hash[:mo] = []
 750             hash[:moinfo] = []
 751           end
 752           while true
 753             #  Loop for each MO
 754             m.clear
 755             ene = nil
 756             spin = nil
 757             sym = nil   #  Not used in Molby
 758             occ = nil
 759             i = 0
 760             while line = getline.call
 761               if line =~ /^ *Sym= *(\w+)/
 762                 sym = $1
 763               elsif line =~ /^ *Ene= *([-+.0-9eE]+)/
 764                 ene = Float($1)
 765               elsif line =~ /^ *Spin= *(\w+)/
 766                 spin = $1
 767               elsif line =~ /^ *Occup= *([-+.0-9eE]+)/
 768                 occ = Float($1)
 769                 if occ > 0.0
 770                   if spin == "Alpha"
 771                     occ_alpha += 1
 772                   else
 773                     occ_beta += 1
 774                   end
 775                 end
 776                 if label
 777                   hash[:moinfo].push([sym, ene, (spin == "Alpha" ? 0 : 1), occ])
 778                 end
 779               elsif line =~ /^ *([0-9]+) +([-+.0-9eE]+)/
 780                 m[i] = Float($2)
 781                 i += 1
 782                 if i >= ncomps
 783                   if spin == "Alpha"
 784                     idx = idx_alpha
 785                     idx_alpha += 1
 786                   else
 787                     idx = idx_beta
 788                     idx_beta += 1
 789                   end
 790                   if label
 791                     hash[:mo].push(m.dup)
 792                   else
 793                     set_mo_coefficients(idx, ene, m)
 794                   end
 795                   break
 796                 end
 797               else
 798                 break
 799               end
 800             end
 801             break if i < ncomps  #  no MO info was found
 802           end
 803           next
 804         end #  end if
 805       end   #  end while
 806     end     #  end catch
 807     if errmsg
 808       message_box("The MOLDEN file was found but not imported. " + errmsg, "Psi4 import info", :ok)
 809       return (label ? nil : false)
 810     end
 811     return (label ? hash : true)
 812   end
 813
 814   #  Import the JANPA log and related molden files
 815   #  Files: inppath.{NAO.molden,CLPO.molden,janpa.log}
 816   #  If inppath.spherical.molden is available, then clear existing mo info
 817   #  and load from it (i.e. use the basis set converted by molden2molden)
 818   def sub_load_janpa_log(inppath)
 819     begin
 820       #  See also: MoleculeGetGaussianComponentInfo() in Molecule.c
 821       m2name_p = {0=>"x", 1=>"y", -1=>"z"}
 822       m2name_d = {0=>"zz-rr", 1=>"xz", -1=>"yz", 2=>"xx-yy", -2=>"xy"}
 823       m2name_f = {0=>"z3-zr2", 1=>"xz2-xr2", -1=>"yz2-yr2", 2=>"x2z-y2z", -2=>"xyz", 3=>"x3-xy2", -3=>"x2y-y3"}
 824       m2name_g = {0=>"z4-z2r2+r4", 1=>"xz3-xzr2", -1=>"yz3-yzr2", 2=>"x2z2-y2z2", -2=>"xyz2-xyr2", 3=>"x3z-xy2z", -3=>"x2yz-y3z", 4=>"x4-x2y2+y4", -4=>"x3y-xy3"}
 825       fp = File.open(inppath + ".janpa.log", "rt") rescue fp = nil
 826       if fp == nil
 827         hide_progress_panel  #  Close if it is open
 828         message_box("Cannot open JANPA log file #{inppath + '.janpa.log'}: " + $!.to_s)
 829         return false
 830       end
 831       print("Importing #{inppath}.janpa.log.\n")
 832       lineno = 0
 833       getline = lambda { lineno += 1; return fp.gets }
 834       h = Hash.new
 835       mfiles = Hash.new
 836       h["software"] = "JANPA"
 837       nao_num = nil  #  Set later
 838       nao_infos = [] #  index=atom_index, value=Hash with key "s", "px", "py" etc.
 839       #  nao_infos[index][key]: array of [nao_num, occupancy], in the reverse order of appearance
 840       while line = getline.call
 841         if line =~ /molden2molden: a conversion tool for MOLDEN/
 842           while line = getline.call
 843             break if line =~ /^All done!/
 844             if line =~ /\.spherical\.molden/
 845               #  The MOs are converted to spherical basis set
 846               #  Clear the existing MO and load *.spherical.molden
 847               sname = inppath + ".spherical.molden"
 848               fps = File.open(sname, "rt") rescue fps = nil
 849               if fps != nil
 850                 print("Importing #{sname}.\n")
 851                 @lineno = 0
 852                 type = get_mo_info(:type)
 853                 alpha = get_mo_info(:alpha)
 854                 beta = get_mo_info(:beta)
 855                 clear_basis_set
 856                 set_mo_info(:type=>type, :alpha=>alpha, :beta=>beta)
 857                 #  mol.@hf_type should be set before calling sub_load_molden
 858                 @hf_type = type
 859                 sub_load_molden(fps)
 860                 fps.close
 861               end
 862             end
 863           end
 864         elsif line =~ /^NAO \#/
 865           h["NAO"] = []  #  [num, anum, occupied_p, group_num, orb_desc, occ, principle*10+orb, shelltype]
 866           # num: NAO index (0-based), anum: atom index (0-based),
 867           # occupied_p: true if marked occupied by JANPA
 868           # group_num: NAO group number (same for same atom, same principle, same orbital type)
 869           # orb_desc: like "s", "px|py|pz", "dxy|dyz|dzx|dx2-y2|dz2-r2", etc.
 870           # occ: occupancy
 871           # principle: principle number (calculated later from group_num)
 872           # orb: 0:s, 1:p, 2:d, 3:f, 4:g
 873           # shelltype: 0: core, 1: valence, 2: rydberg
 874           nao_infos = []  #  nao_infos[atom_index] = {orb_desc=>[[nao0, occ0], [nao1, occ1] ...]}
 875           # orb_desc: same as above
 876           # nao, occ: nao index (0-based) and occupancy
 877           while line = getline.call
 878             break if line !~ /^\s*[1-9]/
 879             num = Integer(line[0, 5]) - 1
 880             name = line[5, 21]
 881             occ = Float(line[26, 11])
 882             #  like A1*: R1*s(0)
 883             #  atom_number, occupied?, group_number, orb_sym, angular_number
 884             name =~ /\s*[A-Z]+([0-9]+)(\*?):\s* R([0-9]+)\*([a-z]+)\(([-0-9]+)\)/
 885             anum = Integer($1) - 1
 886             occupied = $2
 887             group_num = Integer($3)
 888             orb_sym = $4
 889             ang_num = Integer($5)
 890             orb_desc = orb_sym
 891             if orb_desc == "p"
 892               orb_desc += m2name_p[ang_num]
 893             elsif orb_desc == "d"
 894               orb_desc += m2name_d[ang_num]
 895             elsif orb_desc == "f"
 896               orb_desc += m2name_f[ang_num]
 897             elsif orb_desc == "g"
 898               orb_desc += m2name_g[ang_num]
 899             end
 900             h["NAO"].push([num, anum, occupied, group_num, orb_desc, occ, nil, nil])
 901             nao_num = h["NAO"].length
 902             ((nao_infos[anum] ||= Hash.new)[orb_desc] ||= []).unshift([num, occ])
 903           end
 904           nao_num = h["NAO"].length
 905           #  Set principle from nao_infos
 906           val_occ = []  #  val_occ[atom_index][principle*10+orb] = (0: core, 1: valence, 2: ryd)
 907           nao_infos.each_with_index { |value, atom_index|
 908             val_occ[atom_index] ||= []
 909             value.each { |orb_desc, ar|
 910               ar.each_with_index { |v, idx|
 911                 if v[1] > 1.9
 912                   val = 0  #  lp
 913                 elsif v[1] > 0.01
 914                   val = 1  #  valence
 915                 else
 916                   val = 2  #  ryd
 917                 end
 918                 principle = idx + 1
 919                 orb_sym = orb_desc[0]
 920                 orb = 0
 921                 if orb_sym == "p"
 922                   orb = 1
 923                 elsif orb_sym == "d"
 924                   orb = 2
 925                 elsif orb_sym == "f"
 926                   orb = 3
 927                 elsif orb_sym == "g"
 928                   orb = 4
 929                 end
 930                 principle += orb
 931                 po = principle * 10 + orb
 932                 val1 = val_occ[atom_index][po]
 933                 if val1 == nil
 934                   val_occ[atom_index][po] = val
 935                 elsif val1 != val
 936                   val_occ[atom_index][po] = 1  #  core & val, core & ryd or val & ryd
 937                 else
 938                   #  If all orbitals in the shell are "lp", then the shell should be core.
 939                   #  If all orbitals in the shell are "ryd", then the shell should be ryd.
 940                   #  Otherwise, the shell is valence.
 941                 end
 942                 h["NAO"][v[0]][6] = po
 943               }
 944             }
 945           }
 946           nao_num.times { |nao_index|
 947             nao = h["NAO"][nao_index]
 948             nao[7] = val_occ[nao[1]][nao[6]]  #  shell type
 949           }
 950         elsif line =~ /^\s*(C?)LPO\s+D e s c r i p t i o n\s+Occupancy\s+Composition/
 951           if $1 == "C"
 952             key = "CLPO"
 953           else
 954             key = "LPO"
 955           end
 956           h[key] = []  #  [[a1, a2], LP|BD|NB, desc, occ, [h1, h2], lpo_idx, nao_idx]
 957                        # a1, a2: atom indices (0-based), h1, h2: hybrid indices (0-based)
 958                        # lpo_idx: orbital index (0-based)
 959                        # nao_idx: most significant nao (set later)
 960           if key == "CLPO"
 961             h["LHO"] = []  # [[a1, a2], lho_idx, nao_idx]
 962                            # a1, a2: atom indices (0-based)
 963                            # nao_idx: most significant nao (set later)
 964           end
 965           while line = getline.call
 966             break if line =~ /^\s*$/
 967             num = Integer(line[0, 5]) - 1
 968             label1 = line[5, 6].strip
 969             desc = line[11, 30].strip
 970             occ = line[41, 11].strip
 971             comp = line[52, 1000].strip
 972             desc =~ /\s*([-A-Za-z0-9]+)(,\s*(.*$))?/
 973             desc1 = $1
 974             desc2 = ($3 || "")
 975             if desc2 =~ /^(.*)*\(NB\)\s*$/ && label1 == ""
 976               label1 = "(NB)"
 977               desc2 = $1.strip
 978             end
 979             label1.gsub!(/[\(\)]/, "")  #  Strip parens
 980             atoms = desc1.scan(/[A-Za-z]+(\d+)/)   # "C1-H3" -> [["1"], ["3"]]
 981             atoms = atoms.map { |a| Integer(a[0]) - 1 }  # [0, 2]
 982             hybrids_a = comp.scan(/h(\d+)@[A-Za-z]+(\d+)/)  #  "h8@C1...h13@H3" -> "[["8", "1"], ["13", "3"]]
 983             hybrids = []
 984             hybrids_a.each { |a|
 985               i = atoms.find_index(Integer(a[1]) - 1)
 986               if i != nil
 987                 hybrids[i] = Integer(a[0]) - 1
 988               end
 989             } # [7, 12]
 990             #  like ["BD", [0, 2], "Io = 0.2237", occ, [7, 12]]
 991             #  0, 2 are the atom indices (0-based)
 992             #  7, 12 are the number of hybrid orbitals (0-based)
 993             h[key][num] = [atoms, label1, desc2, Float(occ), hybrids, num]
 994             if key == "CLPO"
 995               h["LHO"][hybrids[0]] = [atoms, hybrids[0]]
 996               if hybrids[1]
 997                 h["LHO"][hybrids[1]] = [atoms.reverse, hybrids[1]]
 998               end
 999             end
1000           end
1001           nao_num.times { |idx|
1002             if h[key][idx] == nil
1003               h[key][idx] = [nil, nil, nil, nil, nil, idx]
1004             end
1005           }
1006         elsif line =~ /^ -NAO_Molden_File: (\S*)/
1007           mfiles["NAO"] = $1
1008         elsif line =~ /^ -LHO_Molden_File: (\S*)/
1009           mfiles["LHO"] = $1
1010         elsif line =~ /^ -CLPO_Molden_File: (\S*)/
1011           mfiles["CLPO"] = $1
1012         elsif line =~ /^ -PNAO_Molden_File: (\S*)/
1013           mfiles["PNAO"] = $1
1014         elsif line =~ /^ -AHO_Molden_File: (\S*)/
1015           mfiles["AHO"] = $1
1016         elsif line =~ /^ -LPO_Molden_File: (\S*)/
1017           mfiles["LPO"] = $1
1018         end
1019       end
1020       fp.close
1021
1022       #  Read molden files
1023       keys = []
1024       mfiles.each { |key, value|
1025         fp = Kernel.open(value, "rt") rescue fp = nil
1026         if fp
1027           print("Importing #{value}.\n")
1028           res = sub_load_molden(fp, key)
1029           if res
1030             #  Temporarily assign: this array will be reordered later and converted to LAMatrix
1031             h["AO/#{key}"] = res[:mo]
1032             #h["AO/#{key}"] = LAMatrix.new(res[:mo])
1033           end
1034           fp.close
1035           keys.push(key)
1036         end
1037       }
1038
1039       #  Reorder orbitals
1040       natoms = self.natoms
1041       keys.each { |key|
1042         if key == "NAO"
1043           # [num, anum, occupied_p, group_num, orb_desc, occ, principle*10+orb, shelltype]
1044           # reorder by (1) anum and shelltype (core/val first, ryd later),
1045           # (2) principle*10+orb, (3) num
1046           reorder = []
1047           h["NAO"].each { |nao|
1048             if nao[7] == 2
1049               r = (nao[1] + natoms) * 3 + 2
1050             else
1051               r = nao[1] * 3 + nao[7]
1052             end
1053             r = r * 100 + nao[6]
1054             r = r * nao_num + nao[0]
1055             reorder.push(r)
1056           }
1057           h["NAOnew2old"] = (0...nao_num).to_a.sort { |a, b| reorder[a] <=> reorder[b] }
1058           h["NAOold2new"] = []
1059           h["NAOnew2old"].each_with_index { |i, j|
1060             h["NAOold2new"][i] = j
1061           }
1062         else
1063           #  LPO, CLPO: [[a1, a2], LP|BD|NB, desc, occ, [h1, h2], lpo_idx, nao_idx, shelltype]
1064           #  LHO: [[a1, a2], lho_idx, nao_idx, shelltype]
1065           if h[key]
1066             #  Set the most significant nao
1067             mo = LAMatrix.new(h["AO/NAO"]).inverse * LAMatrix.new(h["AO/#{key}"])
1068             pos = (key == "LHO" ? 2 : 6)  #  position of nao_idx (new index)
1069             nao_num.times { |o_idx|
1070               o = (h[key][o_idx] ||= [])
1071               if o[pos - 1] == nil
1072                 o[pos - 1] = o_idx  #  Original index
1073               end
1074               #  Find the NAO that has the largest contribution to this hybrid orbital
1075               a_max = -1
1076               i_max = nil
1077               nao_num.times { |i|
1078                 a = mo[o_idx, i] ** 2
1079                 if a > a_max
1080                   i_max = i
1081                   a_max = a
1082                 end
1083               }
1084               o[pos] = h["NAOold2new"][i_max]  #  new index
1085               o[pos + 1] = h["NAO"][i_max][7]  #  0:core, 1:val, 2:ryd
1086               if o[0] == nil
1087                 #  Set the atom which has the most significant NAO
1088                 o[0] = [h["NAO"][i_max][1]]
1089               end
1090             }
1091             reorder = []
1092             h[key].each_with_index { |o, o_idx|
1093               #  Reorder by (1) atom and shelltype (core first, val second, ryd later),
1094               #  if LP|BD|NB is present, then val/LP and val/BD come earlier than val/NB
1095               #  (2) second atom if present, and (3) most significant NAO index
1096               if o[pos + 1] == 2
1097                 r = (o[0][0] + natoms) * 4 + 3
1098               else
1099                 r = o[0][0] * 4 + o[pos + 1]
1100                 if o[pos + 1] == 1 && (key == "LPO" || key == "CLPO") && o[1] == "NB"
1101                   r = r + 1
1102                 end
1103               end
1104               r = r * (natoms + 1) + (o[0][1] || -1) + 1  #  Second atom index + 1 (0 if none)
1105               r = r * nao_num + o[pos]
1106               reorder.push(r)
1107             }
1108             h[key + "new2old"] = (0...nao_num).to_a.sort { |a, b| reorder[a] <=> reorder[b] }
1109             h[key + "old2new"] = []
1110             h[key + "new2old"].each_with_index { |i, j|
1111               h[key + "old2new"][i] = j
1112             }
1113           end
1114         end
1115       }
1116       #print("h[\"NAO\"] = " + h["NAO"].inspect + "\n")
1117       #print("h[\"LHO\"] = " + h["LHO"].inspect + "\n")
1118       #print("h[\"LPO\"] = " + h["LPO"].inspect + "\n")
1119       #print("h[\"CLPO\"] = " + h["CLPO"].inspect + "\n")
1120
1121       #  Do reorder and make matrices
1122       keys.each { |key|
1123         if key == "NAO" || key == "PNAO"
1124           old2new = h["NAOold2new"]
1125         elsif key == "LPO"
1126           old2new = h["LPOold2new"]
1127         elsif key == "CLPO"
1128           old2new = h["CLPOold2new"]
1129         elsif key == "AHO" || key == "LHO"
1130           old2new = h["LHOold2new"]
1131         end
1132         if old2new
1133           mo = []
1134           nao_num.times { |i|
1135             mo[old2new[i]] = h["AO/" + key][i]
1136           }
1137           h["AO/" + key] = mo
1138         end
1139         h["AO/" + key] = LAMatrix.new(h["AO/" + key])
1140       }
1141
1142       #  Create labels
1143       keys.each { |key|
1144         old2new = h[key + "old2new"]
1145         if key == "NAO"
1146           h["NAO_L"] = []
1147           nao_num.times { |nao_idx|
1148             # [num, anum, occupied_p, group_num, orb_desc, occ, principle*10+orb, shelltype]
1149             nao = h["NAO"][nao_idx]
1150             aname = self.atoms[nao[1]].name
1151             orb_desc = nao[4]
1152             principle = nao[6] / 10
1153             shell = ["core", "val", "ryd"][nao[7]]
1154             j = (old2new ? old2new[nao_idx] : nao_idx)
1155             h["NAO_L"][j] = "#{aname} (#{principle}#{orb_desc}) (#{shell}, \##{nao_idx + 1})"
1156           }
1157         elsif key == "LHO"
1158           h["LHO_L"] = []
1159           nao_num.times { |lho_idx|
1160             # [[a1, a2], lho_idx, nao_idx, shelltype]
1161             lho = h["LHO"][lho_idx]
1162             aname1 = self.atoms[lho[0][0]].name
1163             aname2 = ("(" + self.atoms[lho[0][1]].name + ")") rescue ""
1164             shell = ["core", "val", "ryd"][lho[3]]
1165             j = (old2new ? old2new[lho_idx] : lho_idx)
1166             h["LHO_L"][j] = "#{aname1}#{aname2} (#{shell}, \##{lho_idx + 1})"
1167           }
1168         elsif key == "LPO" || key == "CLPO"
1169           h[key + "_L"] = []
1170           nao_num.times { |key_idx|
1171             # [[a1, a2], LP|BD|NB, desc, occ, [h1, h2], lpo_idx, nao_idx, shelltype]
1172             o = h[key][key_idx]
1173             aname1 = self.atoms[o[0][0]].name
1174             aname2 = ("-" + self.atoms[o[0][1]].name) rescue ""
1175             type = o[1]
1176             shelltype = o[7]
1177             if shelltype == 2
1178               type = "RY"  #  Rydberg
1179             elsif type == "LP" && shelltype == 0
1180               type = "CR"  #  core
1181             end
1182             j = (old2new ? old2new[key_idx] : key_idx)
1183             h[key + "_L"][j] = "#{aname1}#{aname2} (#{type}, \##{key_idx + 1})"
1184           }
1185         end
1186       }
1187
1188       @nbo = h
1189       if @nbo["AO/NAO"] && @nbo["AO/LHO"] && @nbo["AO/PNAO"]
1190         #  Generate PLHO from PNAO, NAO, LHO
1191         #  This protocol was suggested by the JANPA author in a private commnunication.
1192         begin
1193           nao2lho = @nbo["AO/NAO"].inverse * @nbo["AO/LHO"]
1194           nao2pnao = @nbo["AO/NAO"].inverse * @nbo["AO/PNAO"]
1195           sign = LAMatrix.diagonal((0...nao2pnao.column_size).map { |i| (nao2pnao[i, i] < 0 ? -1 : 1)})
1196           @nbo["AO/PLHO"] = @nbo["AO/PNAO"] * sign * nao2lho
1197         rescue
1198           @nbo["AO/PLHO"] = nil
1199         end
1200       end
1201       if @nbo["AO/AHO"] && @nbo["LHO_L"]
1202         #  Copy labels from LHO
1203         @nbo["AHO_L"] = @nbo["LHO_L"].dup
1204       end
1205       return true
1206     rescue => e
1207       $stderr.write(e.message + "\n")
1208       $stderr.write(e.backtrace.inspect + "\n")
1209     end
1210   end
1211
1212   #  Convert @nbo info to an msbf string (multiline)
1213   #  !:nbo
1214   #  [XXX] (NAO, LHO, etc.)
1215   #  nn (number of components)
1216   #  (nao_labels)*nn, one per line
1217   #  XXX n
1218   #  (%.18g)*nn, at most 6 per line
1219   def nbo2msbfstring
1220     unless @nbo
1221       return ""
1222     end
1223     keys = @nbo.keys.grep(/AO\/\w+/)
1224     s = "!:nbo\n"
1225     keys.each { |k|
1226       k =~ /AO\/(\w+)/
1227       key = $1
1228       labels = @nbo[key + "_L"]
1229       m = @nbo[k]
1230       nc = m.column_size
1231       s += "[#{key}]\n#{nc}\n"
1232       nc.times { |i|
1233         l = (labels ? labels[i] : "")
1234         s += "#{key} #{i+1} #{l}\n"
1235         nc.times { |j|
1236           s += sprintf("%.18g%s", m[i, j], ((j == nc - 1 || j % 6 == 5) ? "\n" : " "))
1237         }
1238       }
1239     }
1240     s += "\n"
1241     s
1242   end
1243
1244   #  Read @nbo info from a multiline string and set to @nbo
1245   def mbsfstring2nbo(s, first_lineno)
1246     h = Hash.new
1247     lineno = first_lineno - 1
1248     errmsg = "Error reading NBO section: "
1249     key = nil
1250     nc = nil
1251     labels = nil
1252     mo_matrix = nil
1253     mo = nil
1254     s.each_line { |ln|
1255       lineno += 1
1256       next if lineno == first_lineno
1257       ln.chomp!
1258       if key == nil
1259         if ln =~ /\[(\w+)\]/
1260           key = $1
1261         else
1262           return errmsg + "keyword (like NBO, NHO, ...) is expected at line #{lineno}"
1263         end
1264         next
1265       end
1266       if nc == nil
1267         nc = ln.to_i
1268         labels = nil
1269         mo_matrix = []
1270         mo = []
1271         next
1272       end
1273       if ln[0, key.length] == key
1274         if mo == nil
1275           return errmsg + "number of components is missing"
1276         end
1277         if mo.length > 0
1278           return errmsg + "too few coefficients at line #{lineno}"
1279         end
1280         ln =~ /^(\w+) +(\d+) *(.*)$/
1281         l = $3
1282         if l != ""
1283           labels ||= []
1284           labels[mo_matrix.length] = l.strip
1285         end
1286         next
1287       end
1288       m = ln.split
1289       begin
1290         m = m.map { |c| Float(c) }
1291       rescue
1292         return errmsg + "cannot convert to number at line #{lineno}"
1293       end
1294       mo += m
1295       if mo.length > nc
1296         return errmsg + "too many coefficients at line #{lineno}"
1297       end
1298       if mo.length == nc
1299         mo_matrix.push(mo)
1300         mo = []
1301         if mo_matrix.length == nc
1302           #  All components are correctly read
1303           if labels
1304             #  Set the label if not specified
1305             nc.times { |i|
1306               if !labels[i] || labels[i] == ""
1307                 labels[i] = "#{key}#{i+1}"
1308               end
1309             }
1310             h["#{key}_L"] = labels
1311           end
1312           h["AO/#{key}"] = LAMatrix.new(mo_matrix)
1313           mo = nil
1314           mo_matrix = nil
1315           labels = nil
1316           key = nil
1317           nc = nil
1318         end
1319       end
1320     }  #  end each_line
1321     @nbo = h
1322     return ""
1323   end
1324
1325   def loadout(filename)
1326   retval = false
1327   fp = open(filename, "rb")
1328   @lineno = 0
1329   begin
1330     while s = fp.gets
1331       @lineno += 1
1332       if s =~ /Gaussian/
1333         retval = sub_load_gaussian_log(fp)
1334         break
1335       elsif s =~ /GAMESS/
1336         retval = sub_load_gamess_log(fp)
1337         break
1338       elsif s =~ /Psi4/
1339         retval = sub_load_psi4_log(fp)
1340         if retval
1341           #  If .molden file exists, then try to read it
1342           namepath = filename.gsub(/\.\w*$/, "")
1343           mname = "#{namepath}.molden"
1344           if File.exists?(mname)
1345             fp2 = open(mname, "rb")
1346             if fp2
1347               flag = sub_load_molden(fp2)
1348               fp2.close
1349               status = (flag ? 0 : -1)
1350             end
1351           end
1352           if File.exists?("#{namepath}.janpa.log")
1353             flag = sub_load_janpa_log(namepath)
1354             status = (flag ? 0 : -1)
1355           end
1356         end
1357         break
1358       end
1359     end
1360   rescue
1361     hide_progress_panel
1362     raise
1363   end
1364         fp.close
1365         return retval
1366   end
1367
1368   alias :loadlog :loadout
1369
1370   def loadxyz(filename)
1371         fp = open(filename, "rb")
1372         n = 0
1373         coords = []
1374         names = []
1375         cell = nil
1376         while 1
1377           line = fp.gets
1378           if line == nil
1379                 fp.close
1380                 break
1381           end
1382           line.chomp
1383           toks = line.split
1384           if coords.length == 0
1385             #  Allow "number of atoms" line or "number of atoms and crystallographic cell parameter"
1386                 #  (Chem3D xyz format)
1387                 next if toks.length == 1
1388                 if toks.length == 7
1389                   cell = toks[1..6].map { |s| Float(s.sub(/\(\d+\)/, "")) }  #  Remove (xx) and convert to float
1390                   next
1391                 end
1392           end
1393           name, x, y, z = line.split
1394           next if z == nil
1395           x = Float(x.sub(/\(\d+\)/, ""))
1396           y = Float(y.sub(/\(\d+\)/, ""))
1397           z = Float(z.sub(/\(\d+\)/, ""))
1398           r = Vector3D[x, y, z]
1399           coords.push(r)
1400           names.push(name)
1401           n += 1
1402         end
1403         celltr = nil
1404         if cell
1405           self.cell = cell
1406           celltr = self.cell_transform
1407         end
1408         names.each_index { |i|
1409           ap = add_atom(names[i])
1410           names[i] =~ /^([A-Za-z]{1,2})/
1411           element = $1.capitalize
1412           ap.element = element
1413           ap.atom_type = element
1414           if celltr
1415             ap.r = celltr * coords[i]
1416           else
1417             ap.r = coords[i]
1418           end
1419         }
1420         guess_bonds
1421         #  Find bonds
1422 #       atoms.each { |ap|
1423 #         j = ap.index
1424 #         (j + 1 ... natoms).each { |k|
1425 #               if calc_bond(j, k) < 1.7
1426 #               #  create_bond(j, k)
1427 #               end
1428 #         }
1429 #       }
1430 #       self.undo_enabled = save_undo_enabled
1431         (n > 0 ? true : false)
1432   end
1433
1434   def savexyz(filename)
1435     open(filename, "wb") { |fp|
1436           fp.printf "%d\n", self.natoms
1437           each_atom { |ap|
1438             fp.printf "%s %.5f %.5f %.5f\n", ap.element, ap.x, ap.y, ap.z
1439           }
1440         }
1441         return true
1442   end
1443
1444   def loadzmat(filename)
1445     self.remove(All)
1446         open(filename, "rb") { |fp|
1447           while (line = fp.gets)
1448             line.chomp!
1449                 a = line.split
1450                 an = Molecule.guess_atomic_number_from_name(a[0])
1451                 elm = Parameter.builtin.elements[an].name
1452                 base1 = a[1].to_i - 1
1453                 base2 = a[3].to_i - 1
1454                 base3 = a[5].to_i - 1
1455                 base1 = nil if base1 < 0
1456                 base2 = nil if base2 < 0
1457                 base3 = nil if base3 < 0
1458                 add_atom(a[0], elm, elm, a[2].to_f, base1, a[4].to_f, base2, a[6].to_f, base3)
1459           end
1460         }
1461         return true
1462   end
1463
1464   def savezmat(filename)
1465   end
1466
1467   def loadcom(filename)
1468 #       save_undo_enabled = self.undo_enabled?
1469 #       self.undo_enabled = false
1470         self.remove(All)
1471         fp = open(filename, "rb")
1472         section = 0
1473         while (line = fp.gets)
1474           line.chomp!
1475           if section == 0
1476             section = 1 if line =~ /^\#/
1477                 next
1478           elsif section == 1 || section == 2
1479             section += 1 if line =~ /^\s*$/
1480                 next
1481           else
1482             #  The first line is skipped (charge and multiplicity)
1483                 while (line = fp.gets)
1484                   line.chomp!
1485                   break if line =~ /^\s*$/
1486                   a = line.split(/\s*[ \t,\/]\s*/)
1487                   r = Vector3D[Float(a[1]), Float(a[2]), Float(a[3])]
1488                   ap = add_atom(a[0])
1489                   a[0] =~ /^([A-Za-z]{1,2})/
1490                   element = $1.capitalize
1491               ap.element = element
1492               ap.atom_type = element
1493               ap.r = r
1494                 end
1495                 break
1496           end
1497         end
1498         fp.close
1499         guess_bonds
1500 #       self.undo_enabled = save_undo_enabled
1501         return true
1502   end
1503
1504   def loadinp(filename)
1505 #       save_undo_enabled = self.undo_enabled?
1506 #       self.undo_enabled = false
1507         self.remove(All)
1508         fp = open(filename, "rb")
1509         section = 0
1510         has_basis = false
1511         while (line = fp.gets)
1512           if line =~ /\A \$BASIS/
1513             has_basis = true
1514                 next
1515           end
1516           next if line !~ /\A \$DATA/
1517           line = fp.gets   #  Title line
1518           line = fp.gets   #  Symmetry line
1519           while (line = fp.gets)
1520 #           puts line
1521             line.chomp!
1522                 break if line =~ /\$END/
1523                 a = line.split
1524                 ap = add_atom(a[0])
1525                 ap.atomic_number = Integer(a[1])
1526                 ap.atom_type = ap.element
1527                 r = Vector3D[Float(a[2]), Float(a[3]), Float(a[4])]
1528                 ap.r = r;
1529                 if !has_basis
1530                   #  Skip until one blank line is detected
1531                   while (line = fp.gets) && line =~ /\S/
1532                   end
1533                 end
1534       end
1535           break
1536         end
1537         fp.close
1538         guess_bonds
1539 #       self.undo_enabled = save_undo_enabled
1540         return true
1541   end
1542
1543   def saveinp(filename)
1544     if natoms == 0
1545       raise MolbyError, "cannot save GAMESS input; the molecule is empty"
1546     end
1547     fp = open(filename, "wb")
1548         now = Time.now.to_s
1549         fp.print <<end_of_header
1550 !  GAMESS input
1551 !  Generated by Molby at #{now}
1552  $CONTRL COORD=UNIQUE EXETYP=RUN ICHARG=0
1553          ICUT=20 INTTYP=HONDO ITOL=30
1554          MAXIT=200 MOLPLT=.T. MPLEVL=0
1555          MULT=1 QMTTOL=1e-08 RUNTYP=OPTIMIZE
1556          SCFTYP=RHF UNITS=ANGS                  $END
1557  $SCF    CONV=1.0E-06 DIRSCF=.T.                $END
1558  $STATPT NSTEP=400 OPTTOL=1.0E-06               $END
1559  $SYSTEM MEMDDI=0 MWORDS=16 TIMLIM=50000        $END
1560  $BASIS  GBASIS=N31 NDFUNC=1 NGAUSS=6           $END
1561  $GUESS  GUESS=HUCKEL                           $END
1562 !
1563  $DATA
1564  #{name}
1565  C1
1566 end_of_header
1567         each_atom { |ap|
1568                 next if ap.atomic_number == 0
1569                 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
1570         }
1571         fp.print " $END\n"
1572         fp.close
1573         return true
1574   end
1575
1576   def savecom(filename)
1577     if natoms == 0
1578       raise MolbyError, "cannot save Gaussian input; the molecule is empty"
1579     end
1580     fp = open(filename, "wb")
1581         base = File.basename(filename, ".*")
1582         fp.print <<end_of_header
1583 %Chk=#{base}.chk
1584 \# PM3 Opt
1585
1586  #{name}; created by Molby at #{Time.now.to_s}
1587
1588  0 1
1589 end_of_header
1590         each_atom { |ap|
1591             next if ap.atomic_number == 0
1592                 fp.printf "%-6s %10.6f %10.6f %10.6f\n", ap.element, ap.r.x, ap.r.y, ap.r.z
1593         }
1594         fp.print "\n"
1595         fp.close
1596         return true
1597   end
1598
1599   alias :loadgjf :loadcom
1600   alias :savegjf :savecom
1601
1602   def loadcif(filename)
1603     mol = self
1604     def getciftoken(fp)
1605           while @tokens.length == 0
1606             line = fp.gets
1607             return nil if !line
1608                 if line[0] == ?;
1609                   s = line
1610                   while line = fp.gets
1611                     break if line[0] == ?;
1612                     s += line
1613                   end
1614                   return s if !line
1615                   s += ";"
1616                   @tokens.push(s)
1617                   line = line[1...line.length]
1618                 end
1619             line.strip!
1620             line.scan(/'[^\']*'|"[^\"]*"|[^#\s]+|#.*/) do |s|
1621               next if s[0] == ?#
1622                   if s =~ /^data_|loop_|global_|save_|stop_/i
1623                     s = "#" + s    #  Label for reserved words
1624                   end
1625                   @tokens.push(s)
1626             end
1627           end
1628           return @tokens.shift
1629         end
1630         def float_strip_rms(str)
1631           str =~ /^(-?)(\d*)(\.(\d*))?(\((\d+)\))?$/
1632           sgn, i, frac, rms = $1, $2, $4, $6
1633           i = i.to_f
1634           if frac
1635                 base = 0.1 ** frac.length
1636                 i = i + frac.to_f * base
1637           else
1638             base = 1.0
1639           end
1640           if rms
1641             rms = rms.to_f * base
1642           else
1643             rms = 0.0
1644           end
1645           if sgn == "-"
1646             i = -i
1647           end
1648           return i, rms
1649         end
1650         def parse_symmetry_operation(str)
1651           if str == "."
1652             sym = nil
1653           elsif (str =~ /(\d+)_(\d)(\d)(\d)/) || (str =~ /(\d+) +(\d)(\d)(\d)/)
1654             sym = [Integer($1) - 1, Integer($2) - 5, Integer($3) - 5, Integer($4) - 5]
1655           elsif (str =~ /^(\d+)$/)
1656             sym = [Integer($1) - 1, 0, 0, 0]
1657           end
1658           if sym && (sym[0] == 0 && sym[1] == 0 && sym[2] == 0 && sym[3] == 0)
1659             sym = nil
1660           end
1661           sym
1662         end
1663         def find_atom_by_name(mol, name)
1664           name = name.delete(" ()")
1665           ap = mol.atoms[name] rescue ap = nil
1666           return ap
1667         end
1668         selfname = self.name
1669         fp = open(filename, "rb")
1670         data_identifier = nil
1671         @tokens = []
1672         count_up = 1
1673         while true
1674           warn_message = ""
1675           verbose = nil
1676           bond_defined = false
1677           special_positions = []
1678           mol.remove(All)
1679           cell = []
1680           cell_trans = cell_trans_inv = Transform.identity
1681           token = getciftoken(fp)
1682           pardigits_re = /\(\d+\)/
1683           calculated_atoms = []
1684           while token != nil
1685             if token =~ /^\#data_/i
1686                   if data_identifier == nil || mol.natoms == 0
1687                     #  First block or no atoms yet
1688             #  Continue processing of this molecule
1689                     data_identifier = token
1690                         token = getciftoken(fp)
1691                         next
1692                   else
1693                     #  Description of another molecule begins here
1694                         data_identifier = token
1695                         break
1696                   end
1697             elsif token =~ /^_cell/
1698                   val = getciftoken(fp)
1699                   if token == "_cell_length_a"
1700                     cell[0], cell[6] = float_strip_rms(val)
1701                   elsif token == "_cell_length_b"
1702                     cell[1], cell[7] = float_strip_rms(val)
1703                   elsif token == "_cell_length_c"
1704                     cell[2], cell[8] = float_strip_rms(val)
1705                   elsif token == "_cell_angle_alpha"
1706                     cell[3], cell[9] = float_strip_rms(val)
1707                   elsif token == "_cell_angle_beta"
1708                     cell[4], cell[10] = float_strip_rms(val)
1709                   elsif token == "_cell_angle_gamma"
1710                     cell[5], cell[11] = float_strip_rms(val)
1711                   end
1712                   if cell.length == 12 && cell.all?
1713                     mol.cell = cell
1714             puts "Unit cell is set to #{mol.cell.inspect}." if verbose
1715                     cell = []
1716                     cell_trans = mol.cell_transform
1717                     cell_trans_inv = cell_trans.inverse
1718                   end
1719                   token = getciftoken(fp)
1720                   next
1721         elsif token.casecmp("#loop_") == 0
1722               labels = []
1723                   while (token = getciftoken(fp)) && token[0] == ?_
1724                     labels.push(token)
1725                   end
1726                   if labels[0] =~ /symmetry_equiv_pos|space_group_symop|atom_site_label|atom_site_aniso_label|geom_bond/
1727                     hlabel = Hash.new(-10000000)
1728                     labels.each_with_index { |lb, i|
1729                           hlabel[lb] = i
1730                     }
1731                     data = []
1732                     n = labels.length
1733                     a = []
1734                     while true
1735                           break if token == nil || token[0] == ?_ || token[0] == ?#
1736                           a.push(token)
1737                           if a.length == n
1738                             data.push(a)
1739                             a = []
1740                           end
1741                           token = getciftoken(fp)
1742                     end
1743                     if labels[0] =~ /^_symmetry_equiv_pos/ || labels[0] =~ /^_space_group_symop/
1744                       data.each { |d|
1745                             symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]] || d[hlabel["_space_group_symop_operation_xyz"]]
1746                             symstr.delete("\"\'")
1747                             exps = symstr.split(/,/)
1748                             sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
1749                             exps.each_with_index { |s, i|
1750                               terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
1751                                   terms.each { |a|
1752                                     #  a[0]: sign, a[2]: numerator, a[4]: denometer
1753                                     if a[4] != nil
1754                                       #  The number part is a[2]/a[4]
1755                                       num = Float(a[2])/Float(a[4])
1756                                     elsif a[2] != nil
1757                                       #  The number part is either integer or a floating point
1758                                       num = Float(a[2])
1759                                     else
1760                                       num = 1.0
1761                                     end
1762                                     num = -num if a[0][0] == ?-
1763                                     xyz = (a[5] || a[6])
1764                                     if xyz == "x" || xyz == "X"
1765                                       sym[i] = num
1766                                     elsif xyz == "y" || xyz == "Y"
1767                                       sym[i + 3] = num
1768                                     elsif xyz == "z" || xyz == "Z"
1769                                       sym[i + 6] = num
1770                                     else
1771                                       sym[9 + i] = num
1772                                     end
1773                                   }
1774                             }
1775                             puts "symmetry operation #{sym.inspect}" if verbose
1776                             mol.add_symmetry(Transform.new(sym))
1777                           }
1778                           puts "#{mol.nsymmetries} symmetry operations are added" if verbose
1779                     elsif labels[0] =~ /^_atom_site_label/
1780                           #  Create atoms
1781                           data.each { |d|
1782                             name = d[hlabel["_atom_site_label"]]
1783                             elem = d[hlabel["_atom_site_type_symbol"]]
1784                             fx = d[hlabel["_atom_site_fract_x"]]
1785                             fy = d[hlabel["_atom_site_fract_y"]]
1786                             fz = d[hlabel["_atom_site_fract_z"]]
1787                             uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
1788                             biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
1789                             occ = d[hlabel["_atom_site_occupancy"]]
1790                             calc = d[hlabel["_atom_site_calc_flag"]]
1791                             name = name.delete(" ()")
1792                             if elem == nil || elem == ""
1793                               if name =~ /[A-Za-z]{1,2}/
1794                                     elem = $&.capitalize
1795                                   else
1796                                     elem = "Du"
1797                                   end
1798                             end
1799                             ap = mol.add_atom(name, elem, elem)
1800                             ap.fract_x, ap.sigma_x = float_strip_rms(fx)
1801                             ap.fract_y, ap.sigma_y = float_strip_rms(fy)
1802                             ap.fract_z, ap.sigma_z = float_strip_rms(fz)
1803                             if biso
1804                               ap.temp_factor, sig = float_strip_rms(biso)
1805                             elsif uiso
1806                               ap.temp_factor, sig = float_strip_rms(uiso)
1807                                   ap.temp_factor *= 78.9568352087149          #  8*pi*pi
1808                             end
1809                             ap.occupancy, sig = float_strip_rms(occ)
1810                             if calc == "c" || calc == "calc"
1811                               calculated_atoms.push(ap.index)
1812                         end
1813                             #  Guess special positions
1814                             (1...mol.nsymmetries).each { |isym|
1815                               sr = ap.fract_r
1816                               sr = (mol.transform_for_symop(isym) * sr) - sr;
1817                                   nx = (sr.x + 0.5).floor
1818                                   ny = (sr.y + 0.5).floor
1819                                   nz = (sr.z + 0.5).floor
1820                                   if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
1821                                     #  [isym, -nx, -ny, -nz] transforms this atom to itself
1822                                     #  The following line is equivalent to:
1823                                     #    if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
1824                                     #    special_positions[ap.index].push(...)
1825                                     (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
1826                                   end
1827                             }
1828                             if verbose && special_positions[ap.index]
1829                               puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
1830                             end
1831                           }
1832                           puts "#{mol.natoms} atoms are created." if verbose
1833                     elsif labels[0] =~ /^_atom_site_aniso_label/
1834                       #  Set anisotropic parameters
1835                           c = 0
1836                           data.each { |d|
1837                             name = d[hlabel["_atom_site_aniso_label"]]
1838                             ap = find_atom_by_name(mol, name)
1839                             next if !ap
1840                             u11 = d[hlabel["_atom_site_aniso_U_11"]]
1841                             if u11
1842                               usig = []
1843                               u11, usig[0] = float_strip_rms(u11)
1844                               u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
1845                               u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
1846                               u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
1847                               u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
1848                               u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
1849                               ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
1850                                   c += 1
1851                             end
1852                           }
1853                           puts "#{c} anisotropic parameters are set." if verbose
1854                     elsif labels[0] =~ /^_geom_bond/
1855                       #  Create bonds
1856                           exbonds = []
1857                           data.each { |d|
1858                             n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
1859                             n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
1860                             sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
1861                             sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
1862                             n1 = find_atom_by_name(mol, n1)
1863                             n2 = find_atom_by_name(mol, n2)
1864                             next if n1 == nil || n2 == nil
1865                         n1 = n1.index
1866                             n2 = n2.index
1867                             sym1 = parse_symmetry_operation(sym1)
1868                             sym2 = parse_symmetry_operation(sym2)
1869                             if sym1 || sym2
1870                               exbonds.push([n1, n2, sym1, sym2])
1871                             else
1872                               mol.create_bond(n1, n2)
1873                             end
1874                             tr1 = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
1875                             tr2 = (sym2 ? mol.transform_for_symop(sym2) : Transform.identity)
1876                             if special_positions[n1]
1877                                   #  Add extra bonds for equivalent positions of n1
1878                                   special_positions[n1].each { |symop|
1879                                     sym2x = mol.symop_for_transform(tr1 * mol.transform_for_symop(symop) * tr1.inverse * tr2)
1880                                     exbonds.push([n1, n2, sym1, sym2x])
1881                                   }
1882                             end
1883                             if special_positions[n2]
1884                                   #  Add extra bonds n2-n1.symop, where symop transforms n2 to self
1885                                   tr = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
1886                                   special_positions[n2].each { |symop|
1887                                     sym1x = mol.symop_for_transform(tr2 * mol.transform_for_symop(symop) * tr2.inverse * tr1)
1888                                     exbonds.push([n2, n1, sym2, sym1x])
1889                                   }
1890                             end
1891                       }
1892               if mol.nbonds > 0
1893                 bond_defined = true
1894               end
1895                           puts "#{mol.nbonds} bonds are created." if verbose
1896                           if calculated_atoms.length > 0
1897                             #  Guess bonds for calculated hydrogen atoms
1898                             n1 = 0
1899                             calculated_atoms.each { |ai|
1900                               if mol.atoms[ai].connects.length == 0
1901                                     as = mol.find_close_atoms(ai)
1902                                     as.each { |aj|
1903                                       mol.create_bond(ai, aj)
1904                                           n1 += 1
1905                                     }
1906                                   end
1907                             }
1908                             puts "#{n1} bonds are guessed." if verbose
1909                           end
1910                           if exbonds.length > 0
1911                             h = Dialog.run("CIF Import: Symmetry Expansion") {
1912                               layout(1,
1913                                     item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
1914                                     item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
1915                                     item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
1916                                     item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
1917                                   )
1918                                   radio_group("atoms_only", "fragment", "ignore")
1919                             }
1920                             if h[:status] == 0 && h["ignore"] == 0
1921                               atoms_only = (h["atoms_only"] != 0)
1922                                   if !atoms_only
1923                                     fragments = []
1924                                     mol.each_fragment { |f| fragments.push(f) }
1925                                   end
1926                                   debug = nil
1927                                   exbonds.each { |ex|
1928                                     if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
1929                                     ex0 = ex.dup
1930                                     (2..3).each { |i|
1931                                       symop = ex[i]
1932                                           if symop == nil
1933                                             ex[i + 2] = ex[i - 2]
1934                                           else
1935                                             if debug; puts "  symop = #{symop.inspect}"; end
1936                                             #  Expand the atom or the fragment including the atom
1937                                             if atoms_only
1938                                                   ig = IntGroup[ex[i - 2]]
1939                                                   idx = 0
1940                                             else
1941                                                   ig = fragments.find { |f| f.include?(ex[i - 2]) }
1942                                                   ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
1943                                             end
1944                                             symop[4] = ex[i - 2]  #  Base atom
1945                                             if debug; puts "  expanding #{ig} by #{symop.inspect}"; end
1946                                             a = mol.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
1947                                             ex[i + 2] = a[idx]   #  Index of the expanded atom
1948                                           end
1949                                     }
1950                                     if ex[4] && ex[5] && ex[4] != ex[5]
1951                                       if debug; puts "  creating bond #{ex[4]} - #{ex[5]}"; end
1952                                       mol.create_bond(ex[4], ex[5])
1953                                     end
1954                                   }
1955                             end
1956                           end
1957                           puts "#{mol.nbonds} bonds are created." if verbose
1958                     end
1959                     next
1960                   else
1961                   #  puts "Loop beginning with #{labels[0]} is skipped"
1962                   end
1963             else
1964               #  Skip this token
1965                   token = getciftoken(fp)
1966             end
1967             #  Skip tokens until next tag or reserved word is detected
1968             while token != nil && token[0] != ?_ && token[0] != ?#
1969                   token = getciftoken(fp)
1970             end
1971             next
1972           end
1973       if !bond_defined
1974                 mol.guess_bonds
1975           end
1976           if token != nil && token == data_identifier
1977             #  Process next molecule: open a new molecule and start adding atom on that
1978                 mol = Molecule.new
1979                 count_up += 1
1980                 (@aux_mols ||= []).push(mol)
1981                 next
1982           end
1983           break
1984         end
1985         fp.close
1986 #       self.undo_enabled = save_undo_enabled
1987         return true
1988   end
1989
1990   def dump(group = nil)
1991     def quote(str)
1992           if str == ""
1993             str = "\"\""
1994           else
1995             str = str.gsub(/%/, "%%")
1996                 str.gsub!(/ /, "%20")
1997                 str.gsub!(/\t/, "%09")
1998           end
1999           str
2000         end
2001         group = atom_group(group ? group : 0...natoms)
2002         s = ""
2003         group.each { |i|
2004           ap = atoms[i]
2005           s += sprintf("%4d %-7s %-4s %-4s %-2s %7.3f %7.3f %7.3f %6.3f [%s]\n",
2006                 ap.index, sprintf("%3s.%d", ap.res_name, ap.res_seq),
2007                 quote(ap.name), quote(ap.atom_type), quote(ap.element),
2008                 ap.r.x, ap.r.y, ap.r.z, ap.charge,
2009                 ap.connects.join(","))
2010         }
2011         print s
2012   end
2013
2014   def from_dump(arg)
2015     if natoms > 0
2016           raise "Molecule must be empty"
2017         end
2018     format = "index residue name atom_type element rx ry rz charge connects"
2019         keys = []
2020         resAtoms = Hash.new
2021         newBonds = []
2022         #  arg can be either a String or an array of String.
2023         #  Iterates for each line in the string or each member of the array.
2024         if arg.is_a?(String)
2025           arg = arg.split("\n")
2026         end
2027     arg.each { |line|
2028           if line =~ /^\#/
2029             format = line[1..-1]
2030                 keys = []
2031           end
2032           if keys.length == 0
2033             keys = format.split(" ").collect { |key| key.to_sym }
2034           end
2035           values = line.chomp.split(" ")
2036           next if values == nil || values.length == 0
2037           ap = create_atom(sprintf("X%03d", natoms))
2038           r = Vector3D[0, 0, 0]
2039           keys.each_index { |i|
2040             break if (value = values[i]) == nil
2041                 if value == "\"\""
2042                   value = ""
2043                 else
2044                   value.gsub(/%09/, "\t")
2045                   value.gsub(/%20/, " ")
2046                   value.gsub(/%%/, "%")
2047                 end
2048                 key = keys[i]
2049                 if key == :residue
2050                   if resAtoms[value] == nil
2051                     resAtoms[value] = []
2052                   end
2053                   resAtoms[value].push(ap.index)
2054                 elsif key == :rx
2055                   r.x = value.to_f
2056                 elsif key == :ry
2057                   r.y = value.to_f
2058                 elsif key == :rz
2059                   r.z = value.to_f
2060                 elsif key == :connects
2061                   value.scan(/\d+/).each { |i|
2062                     i = i.to_i
2063                     if ap.index < i
2064                           newBonds.push(ap.index)
2065                           newBonds.push(i)
2066                         end
2067                   }
2068                 elsif key == :index
2069                   next
2070                 else
2071                   ap.set_attr(key, value)
2072                 end
2073           }
2074           ap.r = r
2075         }
2076         resAtoms.each_key { |key|
2077           assign_residue(atom_group(resAtoms[key]), key)
2078         }
2079         create_bond(*newBonds)
2080         self
2081   end
2082
2083   #  Plug-in for loading mbsf
2084   def loadmbsf_plugin(s, lineno)
2085     ""
2086   end
2087
2088   #  Plug-in for saving mbsf
2089   def savembsf_plugin
2090     ""
2091   end
2092
2093 end