From: Toshi Nagata Date: Sun, 16 Oct 2022 02:04:06 +0000 (+0900) Subject: Start implementing 'additional exponent' for JANPA output X-Git-Url: http://git.osdn.net/view?p=molby%2FMolby.git;a=commitdiff_plain;h=388f1f455139811b1c290d55131fcc1aebf6a5cd Start implementing 'additional exponent' for JANPA output --- diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c index 0dd48b1..4672586 100755 --- a/MolLib/Molecule.c +++ b/MolLib/Molecule.c @@ -1757,7 +1757,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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] == '!') @@ -2710,7 +2710,7 @@ MoleculeLoadShelxFile(Molecule *mp, const char *fname, char **errbuf) /* 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; @@ -2747,6 +2747,7 @@ MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims) 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; diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h index f8ce732..5bc8364 100755 --- a/MolLib/Molecule.h +++ b/MolLib/Molecule.h @@ -219,6 +219,7 @@ typedef struct ShellInfo { 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 */ @@ -410,7 +411,7 @@ Molecule *MoleculeRetain(Molecule *mp); 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); diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c index 6172847..2c09218 100644 --- a/MolLib/Ruby_bind/ruby_bind.c +++ b/MolLib/Ruby_bind/ruby_bind.c @@ -10690,22 +10690,56 @@ s_Molecule_ClearBasisSet(VALUE self) /* * 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) @@ -11882,7 +11916,7 @@ Init_Molby(void) 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); diff --git a/Scripts/gamess.rb b/Scripts/gamess.rb index 1b3dbcd..94c74fe 100755 --- a/Scripts/gamess.rb +++ b/Scripts/gamess.rb @@ -1681,8 +1681,8 @@ class Molecule 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 diff --git a/Scripts/loadsave.rb b/Scripts/loadsave.rb index 728f35a..909d777 100755 --- a/Scripts/loadsave.rb +++ b/Scripts/loadsave.rb @@ -627,6 +627,14 @@ class Molecule 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\]/ @@ -655,61 +663,79 @@ class Molecule 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 @@ -767,8 +793,9 @@ class Molecule 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 @@ -780,27 +807,10 @@ class Molecule # 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 @@ -817,7 +827,30 @@ class Molecule 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]/