+# coding: utf-8
#
# loadsave.rb
#
set_progress_message(mes + "\nReading atomic coordinates...")
if line =~ /ATOMIC/
first_line = true
- if !new_unit
- next # Skip initial atomic coordinates unless loading into an empty molecule
- end
+ # if !new_unit
+ # next # Skip initial atomic coordinates unless loading into an empty molecule
+ # end
line = fp.gets # Skip one line
else
first_line = false
new_unit = false
# create_frame
else
- if search_mode != 1 || nsearch > 1
- # The first frame for geometry search has the same coordinates as input
+ dont_create = false
+ if (search_mode == 1 && nsearch == 1) || first_line
+ # The input coordinate and the first frame for geometry search
+ # can have the same coordinate as the last frame; if this is the case, then
+ # do not create the new frame
+ select_frame(nframes - 1)
+ dont_create = true
+ each_atom { |ap|
+ if (ap.r - coords[ap.index]).length2 > 1e-8
+ dont_create = false
+ break
+ end
+ }
+ end
+ if !dont_create
create_frame([coords]) # Should not be (coords)
end
end
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
(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 "[5D]" or "[9G]".
+ gtos = []
+ spherical = 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 =~ /^\[9G\]/
+ spherical = true
+ # Change the orbital type if necessary
+ gtos.each_with_index { | atom_gtos, atom_index|
+ atom_gtos.each { |gtoline|
+ if gtoline[0] >= 2
+ gtoline[0] = -gtoline[0] # D5/F7/G9
+ end
+ }
+ }
+ 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]
+ 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
alias :savegjf :savecom
def loadcif(filename)
+ mol = self
def getciftoken(fp)
while @tokens.length == 0
line = fp.gets
end
sym
end
- def find_atom_by_name(name)
+ def find_atom_by_name(mol, name)
name = name.delete(" ()")
- ap = self.atoms[name] rescue ap = nil
+ ap = mol.atoms[name] rescue ap = nil
return ap
end
- warn_message = ""
- verbose = nil
- bond_defined = false
- @tokens = []
- special_positions = []
- self.remove(All)
+ selfname = self.name
fp = open(filename, "rb")
- cell = []
- cell_trans = cell_trans_inv = Transform.identity
- token = getciftoken(fp)
- pardigits_re = /\(\d+\)/
- calculated_atoms = []
- while token != nil
- if token =~ /^_cell/
- val = getciftoken(fp)
- if token == "_cell_length_a"
- cell[0], cell[6] = float_strip_rms(val)
- elsif token == "_cell_length_b"
- cell[1], cell[7] = float_strip_rms(val)
- elsif token == "_cell_length_c"
- cell[2], cell[8] = float_strip_rms(val)
- elsif token == "_cell_angle_alpha"
- cell[3], cell[9] = float_strip_rms(val)
- elsif token == "_cell_angle_beta"
- cell[4], cell[10] = float_strip_rms(val)
- elsif token == "_cell_angle_gamma"
- cell[5], cell[11] = float_strip_rms(val)
- end
- if cell.length == 12 && cell.all?
- self.cell = cell
- puts "Unit cell is set to #{cell.inspect}." if verbose
- cell = []
- cell_trans = self.cell_transform
- cell_trans_inv = cell_trans.inverse
- end
- token = getciftoken(fp)
- next
- elsif token.casecmp("#loop_") == 0
- labels = []
- while (token = getciftoken(fp)) && token[0] == ?_
- labels.push(token)
- end
- if labels[0] =~ /symmetry_equiv_pos|atom_site_label|atom_site_aniso_label|geom_bond/
- hlabel = Hash.new(-10000000)
- labels.each_with_index { |lb, i|
- hlabel[lb] = i
- }
- data = []
- n = labels.length
- a = []
- while 1
- break if token == nil || token[0] == ?_ || token[0] == ?#
- a.push(token)
- if a.length == n
- data.push(a)
- a = []
- end
+ data_identifier = nil
+ @tokens = []
+ count_up = 1
+ while true
+ warn_message = ""
+ verbose = nil
+ bond_defined = false
+ special_positions = []
+ mol.remove(All)
+ cell = []
+ cell_trans = cell_trans_inv = Transform.identity
+ token = getciftoken(fp)
+ pardigits_re = /\(\d+\)/
+ calculated_atoms = []
+ while token != nil
+ if token =~ /^\#data_/i
+ if data_identifier == nil || mol.natoms == 0
+ # First block or no atoms yet
+ # Continue processing of this molecule
+ data_identifier = token
token = getciftoken(fp)
+ next
+ else
+ # Description of another molecule begins here
+ data_identifier = token
+ break
+ end
+ elsif token =~ /^_cell/
+ val = getciftoken(fp)
+ if token == "_cell_length_a"
+ cell[0], cell[6] = float_strip_rms(val)
+ elsif token == "_cell_length_b"
+ cell[1], cell[7] = float_strip_rms(val)
+ elsif token == "_cell_length_c"
+ cell[2], cell[8] = float_strip_rms(val)
+ elsif token == "_cell_angle_alpha"
+ cell[3], cell[9] = float_strip_rms(val)
+ elsif token == "_cell_angle_beta"
+ cell[4], cell[10] = float_strip_rms(val)
+ elsif token == "_cell_angle_gamma"
+ cell[5], cell[11] = float_strip_rms(val)
end
- if labels[0] =~ /^_symmetry_equiv_pos/
- data.each { |d|
- symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]]
- symstr.delete("\"\'")
- exps = symstr.split(/,/)
- sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
- exps.each_with_index { |s, i|
- terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
- terms.each { |a|
- # a[0]: sign, a[2]: numerator, a[4]: denometer
- if a[4] != nil
- # The number part is a[2]/a[4]
- num = Float(a[2])/Float(a[4])
- elsif a[2] != nil
- # The number part is either integer or a floating point
- num = Float(a[2])
+ if cell.length == 12 && cell.all?
+ mol.cell = cell
+ puts "Unit cell is set to #{mol.cell.inspect}." if verbose
+ cell = []
+ cell_trans = mol.cell_transform
+ cell_trans_inv = cell_trans.inverse
+ end
+ token = getciftoken(fp)
+ next
+ elsif token.casecmp("#loop_") == 0
+ labels = []
+ while (token = getciftoken(fp)) && token[0] == ?_
+ labels.push(token)
+ end
+ if labels[0] =~ /symmetry_equiv_pos|space_group_symop|atom_site_label|atom_site_aniso_label|geom_bond/
+ hlabel = Hash.new(-10000000)
+ labels.each_with_index { |lb, i|
+ hlabel[lb] = i
+ }
+ data = []
+ n = labels.length
+ a = []
+ while true
+ break if token == nil || token[0] == ?_ || token[0] == ?#
+ a.push(token)
+ if a.length == n
+ data.push(a)
+ a = []
+ end
+ token = getciftoken(fp)
+ end
+ if labels[0] =~ /^_symmetry_equiv_pos/ || labels[0] =~ /^_space_group_symop/
+ data.each { |d|
+ symstr = d[hlabel["_symmetry_equiv_pos_as_xyz"]] || d[hlabel["_space_group_symop_operation_xyz"]]
+ symstr.delete("\"\'")
+ exps = symstr.split(/,/)
+ sym = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
+ exps.each_with_index { |s, i|
+ terms = s.scan(/([-+]?)(([.0-9]+)(\/([0-9]+))?([xXyYzZ])?|([xXyYzZ]))/)
+ terms.each { |a|
+ # a[0]: sign, a[2]: numerator, a[4]: denometer
+ if a[4] != nil
+ # The number part is a[2]/a[4]
+ num = Float(a[2])/Float(a[4])
+ elsif a[2] != nil
+ # The number part is either integer or a floating point
+ num = Float(a[2])
+ else
+ num = 1.0
+ end
+ num = -num if a[0][0] == ?-
+ xyz = (a[5] || a[6])
+ if xyz == "x" || xyz == "X"
+ sym[i] = num
+ elsif xyz == "y" || xyz == "Y"
+ sym[i + 3] = num
+ elsif xyz == "z" || xyz == "Z"
+ sym[i + 6] = num
+ else
+ sym[9 + i] = num
+ end
+ }
+ }
+ puts "symmetry operation #{sym.inspect}" if verbose
+ mol.add_symmetry(Transform.new(sym))
+ }
+ puts "#{mol.nsymmetries} symmetry operations are added" if verbose
+ elsif labels[0] =~ /^_atom_site_label/
+ # Create atoms
+ data.each { |d|
+ name = d[hlabel["_atom_site_label"]]
+ elem = d[hlabel["_atom_site_type_symbol"]]
+ fx = d[hlabel["_atom_site_fract_x"]]
+ fy = d[hlabel["_atom_site_fract_y"]]
+ fz = d[hlabel["_atom_site_fract_z"]]
+ uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
+ biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
+ occ = d[hlabel["_atom_site_occupancy"]]
+ calc = d[hlabel["_atom_site_calc_flag"]]
+ name = name.delete(" ()")
+ if elem == nil || elem == ""
+ if name =~ /[A-Za-z]{1,2}/
+ elem = $&.capitalize
else
- num = 1.0
+ elem = "Du"
end
- num = -num if a[0][0] == ?-
- xyz = (a[5] || a[6])
- if xyz == "x" || xyz == "X"
- sym[i] = num
- elsif xyz == "y" || xyz == "Y"
- sym[i + 3] = num
- elsif xyz == "z" || xyz == "Z"
- sym[i + 6] = num
- else
- sym[9 + i] = num
+ end
+ ap = mol.add_atom(name, elem, elem)
+ ap.fract_x, ap.sigma_x = float_strip_rms(fx)
+ ap.fract_y, ap.sigma_y = float_strip_rms(fy)
+ ap.fract_z, ap.sigma_z = float_strip_rms(fz)
+ if biso
+ ap.temp_factor, sig = float_strip_rms(biso)
+ elsif uiso
+ ap.temp_factor, sig = float_strip_rms(uiso)
+ ap.temp_factor *= 78.9568352087149 # 8*pi*pi
+ end
+ ap.occupancy, sig = float_strip_rms(occ)
+ if calc == "c" || calc == "calc"
+ calculated_atoms.push(ap.index)
+ end
+ # Guess special positions
+ (1...mol.nsymmetries).each { |isym|
+ sr = ap.fract_r
+ sr = (mol.transform_for_symop(isym) * sr) - sr;
+ nx = (sr.x + 0.5).floor
+ ny = (sr.y + 0.5).floor
+ nz = (sr.z + 0.5).floor
+ if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
+ # [isym, -nx, -ny, -nz] transforms this atom to itself
+ # The following line is equivalent to:
+ # if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
+ # special_positions[ap.index].push(...)
+ (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
end
- }
+ }
+ if verbose && special_positions[ap.index]
+ puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
+ end
}
- puts "symmetry operation #{sym.inspect}" if verbose
- add_symmetry(Transform.new(sym))
- }
- puts "#{self.nsymmetries} symmetry operations are added" if verbose
- elsif labels[0] =~ /^_atom_site_label/
- # Create atoms
- data.each { |d|
- name = d[hlabel["_atom_site_label"]]
- elem = d[hlabel["_atom_site_type_symbol"]]
- fx = d[hlabel["_atom_site_fract_x"]]
- fy = d[hlabel["_atom_site_fract_y"]]
- fz = d[hlabel["_atom_site_fract_z"]]
- uiso = d[hlabel["_atom_site_U_iso_or_equiv"]]
- biso = d[hlabel["_atom_site_B_iso_or_equiv"]]
- occ = d[hlabel["_atom_site_occupancy"]]
- calc = d[hlabel["_atom_site_calc_flag"]]
- name = name.delete(" ()")
- if elem == nil || elem == ""
- if name =~ /[A-Za-z]{1,2}/
- elem = $&.capitalize
- else
- elem = "Du"
- end
- end
- ap = self.add_atom(name, elem, elem)
- ap.fract_x, ap.sigma_x = float_strip_rms(fx)
- ap.fract_y, ap.sigma_y = float_strip_rms(fy)
- ap.fract_z, ap.sigma_z = float_strip_rms(fz)
- if biso
- ap.temp_factor, sig = float_strip_rms(biso)
- elsif uiso
- ap.temp_factor, sig = float_strip_rms(uiso)
- ap.temp_factor *= 78.9568352087149 # 8*pi*pi
- end
- ap.occupancy, sig = float_strip_rms(occ)
- if calc == "c" || calc == "calc"
- calculated_atoms.push(ap.index)
- end
- # Guess special positions
- (1...nsymmetries).each { |isym|
- sr = ap.fract_r
- sr = (transform_for_symop(isym) * sr) - sr;
- nx = (sr.x + 0.5).floor
- ny = (sr.y + 0.5).floor
- nz = (sr.z + 0.5).floor
- if (Vector3D[sr.x - nx, sr.y - ny, sr.z - nz].length2 < 1e-6)
- # [isym, -nx, -ny, -nz] transforms this atom to itself
- # The following line is equivalent to:
- # if special_positions[ap.index] == nil; special_positions[ap.index] = []; end;
- # special_positions[ap.index].push(...)
- (special_positions[ap.index] ||= []).push([isym, -nx, -ny, -nz])
- end
+ puts "#{mol.natoms} atoms are created." if verbose
+ elsif labels[0] =~ /^_atom_site_aniso_label/
+ # Set anisotropic parameters
+ c = 0
+ data.each { |d|
+ name = d[hlabel["_atom_site_aniso_label"]]
+ ap = find_atom_by_name(mol, name)
+ next if !ap
+ u11 = d[hlabel["_atom_site_aniso_U_11"]]
+ if u11
+ usig = []
+ u11, usig[0] = float_strip_rms(u11)
+ u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
+ u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
+ u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
+ u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
+ u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
+ ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
+ c += 1
+ end
}
- if verbose && special_positions[ap.index]
- puts "#{name} is on the special position: #{special_positions[ap.index].inspect}"
- end
- }
- puts "#{self.natoms} atoms are created." if verbose
- elsif labels[0] =~ /^_atom_site_aniso_label/
- # Set anisotropic parameters
- c = 0
- data.each { |d|
- name = d[hlabel["_atom_site_aniso_label"]]
- ap = find_atom_by_name(name)
- next if !ap
- u11 = d[hlabel["_atom_site_aniso_U_11"]]
- if u11
- usig = []
- u11, usig[0] = float_strip_rms(u11)
- u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
- u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
- u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
- u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
- u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
- ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
- c += 1
- end
- }
- puts "#{c} anisotropic parameters are set." if verbose
- elsif labels[0] =~ /^_geom_bond/
- # Create bonds
- exbonds = []
- data.each { |d|
- n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
- n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
- sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
- sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
- n1 = find_atom_by_name(n1)
- n2 = find_atom_by_name(n2)
- next if n1 == nil || n2 == nil
- n1 = n1.index
- n2 = n2.index
- sym1 = parse_symmetry_operation(sym1)
- sym2 = parse_symmetry_operation(sym2)
- if sym1 || sym2
- exbonds.push([n1, n2, sym1, sym2])
- else
- self.create_bond(n1, n2)
- end
- tr1 = (sym1 ? transform_for_symop(sym1) : Transform.identity)
- tr2 = (sym2 ? transform_for_symop(sym2) : Transform.identity)
- if special_positions[n1]
- # Add extra bonds for equivalent positions of n1
- special_positions[n1].each { |symop|
- sym2x = symop_for_transform(tr1 * transform_for_symop(symop) * tr1.inverse * tr2)
- exbonds.push([n1, n2, sym1, sym2x])
- }
- end
- if special_positions[n2]
- # Add extra bonds n2-n1.symop, where symop transforms n2 to self
- tr = (sym1 ? transform_for_symop(sym1) : Transform.identity)
- special_positions[n2].each { |symop|
- sym1x = symop_for_transform(tr2 * transform_for_symop(symop) * tr2.inverse * tr1)
- exbonds.push([n2, n1, sym2, sym1x])
- }
- end
- }
- bond_defined = true
- puts "#{self.nbonds} bonds are created." if verbose
- if calculated_atoms.length > 0
- # Guess bonds for calculated hydrogen atoms
- n1 = 0
- calculated_atoms.each { |ai|
- if atoms[ai].connects.length == 0
- as = find_close_atoms(ai)
- as.each { |aj|
- self.create_bond(ai, aj)
- n1 += 1
+ puts "#{c} anisotropic parameters are set." if verbose
+ elsif labels[0] =~ /^_geom_bond/
+ # Create bonds
+ exbonds = []
+ data.each { |d|
+ n1 = d[hlabel["_geom_bond_atom_site_label_1"]]
+ n2 = d[hlabel["_geom_bond_atom_site_label_2"]]
+ sym1 = d[hlabel["_geom_bond_site_symmetry_1"]] || "."
+ sym2 = d[hlabel["_geom_bond_site_symmetry_2"]] || "."
+ n1 = find_atom_by_name(mol, n1)
+ n2 = find_atom_by_name(mol, n2)
+ next if n1 == nil || n2 == nil
+ n1 = n1.index
+ n2 = n2.index
+ sym1 = parse_symmetry_operation(sym1)
+ sym2 = parse_symmetry_operation(sym2)
+ if sym1 || sym2
+ exbonds.push([n1, n2, sym1, sym2])
+ else
+ mol.create_bond(n1, n2)
+ end
+ tr1 = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
+ tr2 = (sym2 ? mol.transform_for_symop(sym2) : Transform.identity)
+ if special_positions[n1]
+ # Add extra bonds for equivalent positions of n1
+ special_positions[n1].each { |symop|
+ sym2x = mol.symop_for_transform(tr1 * mol.transform_for_symop(symop) * tr1.inverse * tr2)
+ exbonds.push([n1, n2, sym1, sym2x])
}
- end
- }
- puts "#{n1} bonds are guessed." if verbose
- end
- if exbonds.length > 0
- h = Dialog.run("CIF Import: Symmetry Expansion") {
- layout(1,
- item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
- item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
- item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
- item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
- )
- radio_group("atoms_only", "fragment", "ignore")
- }
- if h[:status] == 0 && h["ignore"] == 0
- atoms_only = (h["atoms_only"] != 0)
- if !atoms_only
- fragments = []
- self.each_fragment { |f| fragments.push(f) }
- end
- debug = nil
- exbonds.each { |ex|
- if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
- ex0 = ex.dup
- (2..3).each { |i|
- symop = ex[i]
- if symop == nil
- ex[i + 2] = ex[i - 2]
- else
- if debug; puts " symop = #{symop.inspect}"; end
- # Expand the atom or the fragment including the atom
- if atoms_only
- ig = IntGroup[ex[i - 2]]
- idx = 0
+ end
+ if special_positions[n2]
+ # Add extra bonds n2-n1.symop, where symop transforms n2 to self
+ tr = (sym1 ? mol.transform_for_symop(sym1) : Transform.identity)
+ special_positions[n2].each { |symop|
+ sym1x = mol.symop_for_transform(tr2 * mol.transform_for_symop(symop) * tr2.inverse * tr1)
+ exbonds.push([n2, n1, sym2, sym1x])
+ }
+ end
+ }
+ if mol.nbonds > 0
+ bond_defined = true
+ end
+ puts "#{mol.nbonds} bonds are created." if verbose
+ if calculated_atoms.length > 0
+ # Guess bonds for calculated hydrogen atoms
+ n1 = 0
+ calculated_atoms.each { |ai|
+ if mol.atoms[ai].connects.length == 0
+ as = mol.find_close_atoms(ai)
+ as.each { |aj|
+ mol.create_bond(ai, aj)
+ n1 += 1
+ }
+ end
+ }
+ puts "#{n1} bonds are guessed." if verbose
+ end
+ if exbonds.length > 0
+ h = Dialog.run("CIF Import: Symmetry Expansion") {
+ layout(1,
+ item(:text, :title=>"There are bonds including symmetry related atoms.\nWhat do you want to do?"),
+ item(:radio, :title=>"Expand only atoms that are included in those extra bonds.", :tag=>"atoms_only"),
+ item(:radio, :title=>"Expand fragments having atoms included in the extra bonds.", :tag=>"fragment", :value=>1),
+ item(:radio, :title=>"Ignore these extra bonds.", :tag=>"ignore")
+ )
+ radio_group("atoms_only", "fragment", "ignore")
+ }
+ if h[:status] == 0 && h["ignore"] == 0
+ atoms_only = (h["atoms_only"] != 0)
+ if !atoms_only
+ fragments = []
+ mol.each_fragment { |f| fragments.push(f) }
+ end
+ debug = nil
+ exbonds.each { |ex|
+ if debug; puts "extra bond #{ex[0]}(#{ex[2].inspect}) - #{ex[1]}(#{ex[3].inspect})"; end
+ ex0 = ex.dup
+ (2..3).each { |i|
+ symop = ex[i]
+ if symop == nil
+ ex[i + 2] = ex[i - 2]
else
- ig = fragments.find { |f| f.include?(ex[i - 2]) }
- ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
+ if debug; puts " symop = #{symop.inspect}"; end
+ # Expand the atom or the fragment including the atom
+ if atoms_only
+ ig = IntGroup[ex[i - 2]]
+ idx = 0
+ else
+ ig = fragments.find { |f| f.include?(ex[i - 2]) }
+ ig.each_with_index { |n, ii| if n == ex[i - 2]; idx = ii; break; end }
+ end
+ symop[4] = ex[i - 2] # Base atom
+ if debug; puts " expanding #{ig} by #{symop.inspect}"; end
+ a = mol.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
+ ex[i + 2] = a[idx] # Index of the expanded atom
end
- symop[4] = ex[i - 2] # Base atom
- if debug; puts " expanding #{ig} by #{symop.inspect}"; end
- a = self.expand_by_symmetry(ig, symop[0], symop[1], symop[2], symop[3])
- ex[i + 2] = a[idx] # Index of the expanded atom
- end
+ }
+ if ex[4] && ex[5] && ex[4] != ex[5]
+ if debug; puts " creating bond #{ex[4]} - #{ex[5]}"; end
+ mol.create_bond(ex[4], ex[5])
+ end
}
- if ex[4] && ex[5] && ex[4] != ex[5]
- if debug; puts " creating bond #{ex[4]} - #{ex[5]}"; end
- self.create_bond(ex[4], ex[5])
- end
- }
+ end
end
- end
- puts "#{self.nbonds} bonds are created." if verbose
+ puts "#{mol.nbonds} bonds are created." if verbose
+ end
+ next
+ else
+ # puts "Loop beginning with #{labels[0]} is skipped"
end
- next
- else
- # puts "Loop beginning with #{labels[0]} is skipped"
- end
- else
- # Skip this token
- token = getciftoken(fp)
+ else
+ # Skip this token
+ token = getciftoken(fp)
+ end
+ # Skip tokens until next tag or reserved word is detected
+ while token != nil && token[0] != ?_ && token[0] != ?#
+ token = getciftoken(fp)
+ end
+ next
+ end
+ if !bond_defined
+ mol.guess_bonds
end
- # Skip tokens until next tag or reserved word is detected
- while token != nil && token[0] != ?_ && token[0] != ?#
- token = getciftoken(fp)
+ if token != nil && token == data_identifier
+ # Process next molecule: open a new molecule and start adding atom on that
+ mol = Molecule.new
+ count_up += 1
+ (@aux_mols ||= []).push(mol)
+ next
end
- next
+ break
end
fp.close
- if !bond_defined
- self.guess_bonds
- end
# self.undo_enabled = save_undo_enabled
return true
end
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 = []
self
end
+ # Plug-in for loading mbsf
+ def loadmbsf_plugin(s, lineno)
+ ""
+ end
+
+ # Plug-in for saving mbsf
+ def savembsf_plugin
+ ""
+ end
+
end