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