s_append_asprintf(errbuf, "line %d: the atom index (%d) is out of range", lineNumber, ibuf[1]);
goto err_exit;
}
- MoleculeAddGaussianOrbitalShell(mp, ibuf[1], ibuf[2], ibuf[0]);
+ MoleculeAddGaussianOrbitalShell(mp, ibuf[1], ibuf[2], ibuf[0], 0);
i = ibuf[0];
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
/* Add one gaussian orbital shell information (not undoable) */
int
-MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims)
+MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims, Int add_exp)
{
BasisSet *bset;
ShellInfo *shellp;
shellp->m_idx = 0;
shellp->p_idx = 0;
}
+ shellp->add_exp = add_exp;
/* Update the number of components (if not yet determined) */
if (bset->ncomps < shellp->m_idx + shellp->ncomp)
bset->ncomps = shellp->m_idx + shellp->ncomp;
signed char sym; /* Symmetry of the basis; S, P, ... */
signed char ncomp; /* Number of components (S: 1, P: 3, SP: 4, etc.) */
signed char nprim; /* Number of primitives for this shell */
+ signed char add_exp; /* Additional exponent (for JANPA-Molden only) */
Int p_idx; /* Index to the PrimInfo (exponent/coefficient) table */
Int cn_idx; /* Index to the normalized (cached) contraction coefficient table */
Int a_idx; /* Index to the atom which this primitive belongs to */
void MoleculeRelease(Molecule *mp);
void MoleculeExchange(Molecule *mp1, Molecule *mp2);
-int MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims);
+int MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims, Int add_exp);
int MoleculeAddGaussianPrimitiveCoefficients(Molecule *mol, Double exponent, Double contraction, Double contraction_sp);
int MoleculeGetGaussianComponentInfo(Molecule *mol, Int comp_idx, Int *outAtomIdx, char *outLabel, Int *outShellIdx);
int MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Double *coeffs);
/*
* call-seq:
- * add_gaussian_orbital_shell(atom_index, sym, no_of_primitives)
+ * add_gaussian_orbital_shell(atom_index, sym, no_of_primitives[, additional_exponent])
*
* To be used internally. Add a gaussian orbital shell with the atom index, symmetry code,
- * and the number of primitives. Symmetry code: 0, S-type; 1, P-type; -1, SP-type; 2, D-type;
- * -2, D5-type.
+ * and the number of primitives.
+ * Additional exponent is for JANPA only; implements an additinal r^N component that
+ * appears in cartesian->spherical conversion.
+ * Symmetry code: 0, S-type; 1, P-type; -1, SP-type; 2, D-type; -2, D5-type;
+ * 3, F-type; -3, F7-type; 4, G-type; -4, G9-type.
+ * Or: "s", S-type; "p", P-type; "sp", SP-type; "d", D-type; "d5", D5-type;
+ * "f", F-type; "f7", F7-type; "g", G-type; "g9", G9-type
*/
static VALUE
-s_Molecule_AddGaussianOrbitalShell(VALUE self, VALUE aval, VALUE symval, VALUE npval)
+s_Molecule_AddGaussianOrbitalShell(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- int sym, nprims, a_idx, n;
- Data_Get_Struct(self, Molecule, mol);
+ int sym, nprims, a_idx, n, add_exp;
+ VALUE aval, symval, npval, addval;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "31", &aval, &symval, &npval, &addval);
+ if (rb_obj_is_kind_of(symval, rb_cString)) {
+ const char *p = StringValuePtr(symval);
+ if (strcasecmp(p, "s") == 0)
+ sym = 0;
+ else if (strcasecmp(p, "p") == 0)
+ sym = 1;
+ else if (strcasecmp(p, "sp") == 0)
+ sym = -1;
+ else if (strcasecmp(p, "d") == 0)
+ sym = 2;
+ else if (strcasecmp(p, "d5") == 0)
+ sym = -2;
+ else if (strcasecmp(p, "f") == 0)
+ sym = 3;
+ else if (strcasecmp(p, "f7") == 0)
+ sym = -3;
+ else if (strcasecmp(p, "g") == 0)
+ sym = 4;
+ else if (strcasecmp(p, "g9") == 0)
+ sym = -4;
+ else
+ rb_raise(rb_eArgError, "Unknown orbital type '%s'", p);
+ } else {
+ sym = NUM2INT(rb_Integer(symval));
+ }
a_idx = NUM2INT(rb_Integer(aval));
- sym = NUM2INT(rb_Integer(symval));
nprims = NUM2INT(rb_Integer(npval));
- n = MoleculeAddGaussianOrbitalShell(mol, a_idx, sym, nprims);
+ if (addval != Qnil)
+ add_exp = NUM2INT(rb_Integer(addval));
+ else add_exp = 0;
+ n = MoleculeAddGaussianOrbitalShell(mol, a_idx, sym, nprims, add_exp);
if (n == -1)
rb_raise(rb_eMolbyError, "Molecule is emptry");
else if (n == -2)
rb_define_method(rb_cMolecule, "nelpots", s_Molecule_NElpots, 0);
rb_define_method(rb_cMolecule, "elpot", s_Molecule_Elpot, 1);
rb_define_method(rb_cMolecule, "clear_basis_set", s_Molecule_ClearBasisSet, 0);
- rb_define_method(rb_cMolecule, "add_gaussian_orbital_shell", s_Molecule_AddGaussianOrbitalShell, 3);
+ rb_define_method(rb_cMolecule, "add_gaussian_orbital_shell", s_Molecule_AddGaussianOrbitalShell, -1);
rb_define_method(rb_cMolecule, "add_gaussian_primitive_coefficients", s_Molecule_AddGaussianPrimitiveCoefficients, 3);
rb_define_method(rb_cMolecule, "get_gaussian_shell_info", s_Molecule_GetGaussianShellInfo, 1);
rb_define_method(rb_cMolecule, "get_gaussian_primitive_coefficients", s_Molecule_GetGaussianPrimitiveCoefficients, 1);
if flag
if mol
# import JANPA log and molden output
- # Files: inppath.janpa.log, inppath.{NAO,PNAO,AHO,LHO,LPO,CLPO}.molden
- mol.sub_load_janpa_log(inppath, spherical)
+ # Files: inppath.janpa.log, inppath.{NAO,PNAO,AHO,LHO,LPO,CLPO,spherical}.molden
+ mol.sub_load_janpa_log(inppath)
end
hide_progress_panel
else
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\]/
elsif line =~ /^\[GTO\]/
shell = 0
atom_index = 0
- if label
- gtos = []
- hash[:gto] = []
- end
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
- case a[0]
+ a[0] =~ /^([a-z]+)([0-9]+)?$/
+ symcode = $1
+ add_exp = ($2 == nil ? 0 : $2.to_i)
+ case symcode
when "s"
- sym = 0; n = 1
+ sym = 0
when "p"
- sym = 1; n = 3
+ sym = 1
when "d"
- sym = 2; n = 6 # TODO: handle both spherical and cartesian
+ sym = 2
when "f"
- sym = 3; n = 10
+ sym = 3
when "g"
- sym = 4; n = 15
+ sym = 4
else
raise MolbyError, "Unknown gaussian shell type '#{a[0]}' at line #{@lineno} in MOLDEN file"
end
nprimitives = Integer(a[1])
- if label
- gtoline = [sym, []]
- gtos.push(gtoline)
- end
+ gtoline = [sym, [], add_exp]
+ atom_gtos.push(gtoline)
nprimitives.times { |i|
line = getline.call # exponent, contraction
b = line.split(' ')
- if label
- gtoline[1].push(Float(b[0]), Float(b[1]))
- end
- add_gaussian_primitive_coefficients(Float(b[0]), Float(b[1]), 0.0)
+ gtoline[1].push(Float(b[0]), Float(b[1]))
}
# end of one shell
- if label == nil
- add_gaussian_orbital_shell(atom_index, sym, nprimitives)
- end
shell += 1
- ncomps += n
end
# end of one atom
atom_index += 1
- if label
- hash[:gto].push(gtos)
- end
+ 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
end
break if i < ncomps # no MO info was found
end
+ # TODO: reorder D, F, G coefficients for Molby order
next
- end
+ end # end if
end # end while
end # end catch
if errmsg
# Import the JANPA log and related molden files
# Files: inppath.{NAO.molden,CLPO.molden,janpa.log}
- # If mo_path is not nil, then clear existing mo info and load from mo_path
- def sub_load_janpa_log(inppath, mo_path = nil)
+ # 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
- if mo_path
- fp = File.open(mo_path, "rt") rescue fp = nil
- if fp == nil
- hide_progress_panel # Close if it is open
- message_box("Cannot open MO molden file #{mo_path}: " + $!.to_s)
- return false
- end
- @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(fp)
- fp.close
- end
fp = File.open(inppath + ".janpa.log", "rt") rescue fp = nil
if fp == nil
hide_progress_panel # Close if it is open
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 =~ /^NAO \#/
+ 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]/
OSDN Git Service
