(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
calculated_atoms = []
while token != nil
if token =~ /^\#data_/i
- if token.casecmp("#data_global") == 0
- token = getciftoken(fp)
- next
- elsif data_identifier == nil
- # Continue processing of this molecule
+ 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
end
if cell.length == 12 && cell.all?
mol.cell = cell
- puts "Unit cell is set to #{cell.inspect}." if verbose
+ puts "Unit cell is set to #{mol.cell.inspect}." if verbose
cell = []
- cell_trans = self.cell_transform
+ cell_trans = mol.cell_transform
cell_trans_inv = cell_trans.inverse
end
token = getciftoken(fp)
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/
+ 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
end
token = getciftoken(fp)
end
- if labels[0] =~ /^_symmetry_equiv_pos/
+ if labels[0] =~ /^_symmetry_equiv_pos/ || labels[0] =~ /^_space_group_symop/
data.each { |d|
- symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]]
+ 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]
}
end
}
- bond_defined = true
+ 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
OSDN Git Service
