Molecular Modeling Software
リビジョン | 388f1f455139811b1c290d55131fcc1aebf6a5cd (tree) |
---|---|
日時 | 2022-10-16 11:04:06 |
作者 | Toshi Nagata <alchemist.2005@nift...> |
コミッター | Toshi Nagata |
Start implementing 'additional exponent' for JANPA output
@@ -1757,7 +1757,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) | ||
1757 | 1757 | s_append_asprintf(errbuf, "line %d: the atom index (%d) is out of range", lineNumber, ibuf[1]); |
1758 | 1758 | goto err_exit; |
1759 | 1759 | } |
1760 | - MoleculeAddGaussianOrbitalShell(mp, ibuf[1], ibuf[2], ibuf[0]); | |
1760 | + MoleculeAddGaussianOrbitalShell(mp, ibuf[1], ibuf[2], ibuf[0], 0); | |
1761 | 1761 | i = ibuf[0]; |
1762 | 1762 | while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) { |
1763 | 1763 | if (buf[0] == '!') |
@@ -2710,7 +2710,7 @@ MoleculeLoadShelxFile(Molecule *mp, const char *fname, char **errbuf) | ||
2710 | 2710 | |
2711 | 2711 | /* Add one gaussian orbital shell information (not undoable) */ |
2712 | 2712 | int |
2713 | -MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims) | |
2713 | +MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims, Int add_exp) | |
2714 | 2714 | { |
2715 | 2715 | BasisSet *bset; |
2716 | 2716 | ShellInfo *shellp; |
@@ -2747,6 +2747,7 @@ MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims) | ||
2747 | 2747 | shellp->m_idx = 0; |
2748 | 2748 | shellp->p_idx = 0; |
2749 | 2749 | } |
2750 | + shellp->add_exp = add_exp; | |
2750 | 2751 | /* Update the number of components (if not yet determined) */ |
2751 | 2752 | if (bset->ncomps < shellp->m_idx + shellp->ncomp) |
2752 | 2753 | bset->ncomps = shellp->m_idx + shellp->ncomp; |
@@ -219,6 +219,7 @@ typedef struct ShellInfo { | ||
219 | 219 | signed char sym; /* Symmetry of the basis; S, P, ... */ |
220 | 220 | signed char ncomp; /* Number of components (S: 1, P: 3, SP: 4, etc.) */ |
221 | 221 | signed char nprim; /* Number of primitives for this shell */ |
222 | + signed char add_exp; /* Additional exponent (for JANPA-Molden only) */ | |
222 | 223 | Int p_idx; /* Index to the PrimInfo (exponent/coefficient) table */ |
223 | 224 | Int cn_idx; /* Index to the normalized (cached) contraction coefficient table */ |
224 | 225 | Int a_idx; /* Index to the atom which this primitive belongs to */ |
@@ -410,7 +411,7 @@ Molecule *MoleculeRetain(Molecule *mp); | ||
410 | 411 | void MoleculeRelease(Molecule *mp); |
411 | 412 | void MoleculeExchange(Molecule *mp1, Molecule *mp2); |
412 | 413 | |
413 | -int MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims); | |
414 | +int MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims, Int add_exp); | |
414 | 415 | int MoleculeAddGaussianPrimitiveCoefficients(Molecule *mol, Double exponent, Double contraction, Double contraction_sp); |
415 | 416 | int MoleculeGetGaussianComponentInfo(Molecule *mol, Int comp_idx, Int *outAtomIdx, char *outLabel, Int *outShellIdx); |
416 | 417 | int MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Double *coeffs); |
@@ -10690,22 +10690,56 @@ s_Molecule_ClearBasisSet(VALUE self) | ||
10690 | 10690 | |
10691 | 10691 | /* |
10692 | 10692 | * call-seq: |
10693 | - * add_gaussian_orbital_shell(atom_index, sym, no_of_primitives) | |
10693 | + * add_gaussian_orbital_shell(atom_index, sym, no_of_primitives[, additional_exponent]) | |
10694 | 10694 | * |
10695 | 10695 | * To be used internally. Add a gaussian orbital shell with the atom index, symmetry code, |
10696 | - * and the number of primitives. Symmetry code: 0, S-type; 1, P-type; -1, SP-type; 2, D-type; | |
10697 | - * -2, D5-type. | |
10696 | + * and the number of primitives. | |
10697 | + * Additional exponent is for JANPA only; implements an additinal r^N component that | |
10698 | + * appears in cartesian->spherical conversion. | |
10699 | + * Symmetry code: 0, S-type; 1, P-type; -1, SP-type; 2, D-type; -2, D5-type; | |
10700 | + * 3, F-type; -3, F7-type; 4, G-type; -4, G9-type. | |
10701 | + * Or: "s", S-type; "p", P-type; "sp", SP-type; "d", D-type; "d5", D5-type; | |
10702 | + * "f", F-type; "f7", F7-type; "g", G-type; "g9", G9-type | |
10698 | 10703 | */ |
10699 | 10704 | static VALUE |
10700 | -s_Molecule_AddGaussianOrbitalShell(VALUE self, VALUE aval, VALUE symval, VALUE npval) | |
10705 | +s_Molecule_AddGaussianOrbitalShell(int argc, VALUE *argv, VALUE self) | |
10701 | 10706 | { |
10702 | 10707 | Molecule *mol; |
10703 | - int sym, nprims, a_idx, n; | |
10704 | - Data_Get_Struct(self, Molecule, mol); | |
10708 | + int sym, nprims, a_idx, n, add_exp; | |
10709 | + VALUE aval, symval, npval, addval; | |
10710 | + Data_Get_Struct(self, Molecule, mol); | |
10711 | + rb_scan_args(argc, argv, "31", &aval, &symval, &npval, &addval); | |
10712 | + if (rb_obj_is_kind_of(symval, rb_cString)) { | |
10713 | + const char *p = StringValuePtr(symval); | |
10714 | + if (strcasecmp(p, "s") == 0) | |
10715 | + sym = 0; | |
10716 | + else if (strcasecmp(p, "p") == 0) | |
10717 | + sym = 1; | |
10718 | + else if (strcasecmp(p, "sp") == 0) | |
10719 | + sym = -1; | |
10720 | + else if (strcasecmp(p, "d") == 0) | |
10721 | + sym = 2; | |
10722 | + else if (strcasecmp(p, "d5") == 0) | |
10723 | + sym = -2; | |
10724 | + else if (strcasecmp(p, "f") == 0) | |
10725 | + sym = 3; | |
10726 | + else if (strcasecmp(p, "f7") == 0) | |
10727 | + sym = -3; | |
10728 | + else if (strcasecmp(p, "g") == 0) | |
10729 | + sym = 4; | |
10730 | + else if (strcasecmp(p, "g9") == 0) | |
10731 | + sym = -4; | |
10732 | + else | |
10733 | + rb_raise(rb_eArgError, "Unknown orbital type '%s'", p); | |
10734 | + } else { | |
10735 | + sym = NUM2INT(rb_Integer(symval)); | |
10736 | + } | |
10705 | 10737 | a_idx = NUM2INT(rb_Integer(aval)); |
10706 | - sym = NUM2INT(rb_Integer(symval)); | |
10707 | 10738 | nprims = NUM2INT(rb_Integer(npval)); |
10708 | - n = MoleculeAddGaussianOrbitalShell(mol, a_idx, sym, nprims); | |
10739 | + if (addval != Qnil) | |
10740 | + add_exp = NUM2INT(rb_Integer(addval)); | |
10741 | + else add_exp = 0; | |
10742 | + n = MoleculeAddGaussianOrbitalShell(mol, a_idx, sym, nprims, add_exp); | |
10709 | 10743 | if (n == -1) |
10710 | 10744 | rb_raise(rb_eMolbyError, "Molecule is emptry"); |
10711 | 10745 | else if (n == -2) |
@@ -11882,7 +11916,7 @@ Init_Molby(void) | ||
11882 | 11916 | rb_define_method(rb_cMolecule, "nelpots", s_Molecule_NElpots, 0); |
11883 | 11917 | rb_define_method(rb_cMolecule, "elpot", s_Molecule_Elpot, 1); |
11884 | 11918 | rb_define_method(rb_cMolecule, "clear_basis_set", s_Molecule_ClearBasisSet, 0); |
11885 | - rb_define_method(rb_cMolecule, "add_gaussian_orbital_shell", s_Molecule_AddGaussianOrbitalShell, 3); | |
11919 | + rb_define_method(rb_cMolecule, "add_gaussian_orbital_shell", s_Molecule_AddGaussianOrbitalShell, -1); | |
11886 | 11920 | rb_define_method(rb_cMolecule, "add_gaussian_primitive_coefficients", s_Molecule_AddGaussianPrimitiveCoefficients, 3); |
11887 | 11921 | rb_define_method(rb_cMolecule, "get_gaussian_shell_info", s_Molecule_GetGaussianShellInfo, 1); |
11888 | 11922 | rb_define_method(rb_cMolecule, "get_gaussian_primitive_coefficients", s_Molecule_GetGaussianPrimitiveCoefficients, 1); |
@@ -1681,8 +1681,8 @@ class Molecule | ||
1681 | 1681 | if flag |
1682 | 1682 | if mol |
1683 | 1683 | # import JANPA log and molden output |
1684 | - # Files: inppath.janpa.log, inppath.{NAO,PNAO,AHO,LHO,LPO,CLPO}.molden | |
1685 | - mol.sub_load_janpa_log(inppath, spherical) | |
1684 | + # Files: inppath.janpa.log, inppath.{NAO,PNAO,AHO,LHO,LPO,CLPO,spherical}.molden | |
1685 | + mol.sub_load_janpa_log(inppath) | |
1686 | 1686 | end |
1687 | 1687 | hide_progress_panel |
1688 | 1688 | else |
@@ -627,6 +627,14 @@ class Molecule | ||
627 | 627 | if label |
628 | 628 | hash = Hash.new |
629 | 629 | end |
630 | + # The GTOs (orbital type, contractions and exponents) are stored in gtos[] | |
631 | + # and set just before first [MO] is processed. | |
632 | + # This is because we do not know whether the orbital type is cartesian or spherical | |
633 | + # until we see "[5D]" or "[9G]". | |
634 | + gtos = [] | |
635 | + spherical = false | |
636 | + # Number of components for each orbital type | |
637 | + ncomp_hash = { 0=>1, 1=>3, -1=>4, 2=>6, -2=>5, 3=>10, -3=>7, 4=>15, -4=>9 } | |
630 | 638 | catch :ignore do |
631 | 639 | while line = getline.call |
632 | 640 | if line =~ /^\[Atoms\]/ |
@@ -655,61 +663,79 @@ class Molecule | ||
655 | 663 | elsif line =~ /^\[GTO\]/ |
656 | 664 | shell = 0 |
657 | 665 | atom_index = 0 |
658 | - if label | |
659 | - gtos = [] | |
660 | - hash[:gto] = [] | |
661 | - end | |
662 | 666 | while line = getline.call |
663 | 667 | # index, 0? |
664 | 668 | a = line.split(' ') |
665 | 669 | break if a.length != 2 |
670 | + atom_gtos = [] # [[sym1, [e11, c11, e12, c12, ...], add_exp1], [sym2, [e21, c22, ...], add_exp2], ...] | |
666 | 671 | # loop for shells |
667 | 672 | while line = getline.call |
668 | 673 | # type, no_of_primitives, 1.00? |
669 | 674 | a = line.split(' ') |
670 | 675 | break if a.length != 3 # Terminated by a blank line |
671 | - case a[0] | |
676 | + a[0] =~ /^([a-z]+)([0-9]+)?$/ | |
677 | + symcode = $1 | |
678 | + add_exp = ($2 == nil ? 0 : $2.to_i) | |
679 | + case symcode | |
672 | 680 | when "s" |
673 | - sym = 0; n = 1 | |
681 | + sym = 0 | |
674 | 682 | when "p" |
675 | - sym = 1; n = 3 | |
683 | + sym = 1 | |
676 | 684 | when "d" |
677 | - sym = 2; n = 6 # TODO: handle both spherical and cartesian | |
685 | + sym = 2 | |
678 | 686 | when "f" |
679 | - sym = 3; n = 10 | |
687 | + sym = 3 | |
680 | 688 | when "g" |
681 | - sym = 4; n = 15 | |
689 | + sym = 4 | |
682 | 690 | else |
683 | 691 | raise MolbyError, "Unknown gaussian shell type '#{a[0]}' at line #{@lineno} in MOLDEN file" |
684 | 692 | end |
685 | 693 | nprimitives = Integer(a[1]) |
686 | - if label | |
687 | - gtoline = [sym, []] | |
688 | - gtos.push(gtoline) | |
689 | - end | |
694 | + gtoline = [sym, [], add_exp] | |
695 | + atom_gtos.push(gtoline) | |
690 | 696 | nprimitives.times { |i| |
691 | 697 | line = getline.call # exponent, contraction |
692 | 698 | b = line.split(' ') |
693 | - if label | |
694 | - gtoline[1].push(Float(b[0]), Float(b[1])) | |
695 | - end | |
696 | - add_gaussian_primitive_coefficients(Float(b[0]), Float(b[1]), 0.0) | |
699 | + gtoline[1].push(Float(b[0]), Float(b[1])) | |
697 | 700 | } |
698 | 701 | # end of one shell |
699 | - if label == nil | |
700 | - add_gaussian_orbital_shell(atom_index, sym, nprimitives) | |
701 | - end | |
702 | 702 | shell += 1 |
703 | - ncomps += n | |
704 | 703 | end |
705 | 704 | # end of one atom |
706 | 705 | atom_index += 1 |
707 | - if label | |
708 | - hash[:gto].push(gtos) | |
709 | - end | |
706 | + gtos.push(atom_gtos) | |
707 | + end | |
708 | + if label | |
709 | + hash[:gto] = gtos | |
710 | 710 | end |
711 | 711 | redo # The next line will be the beginning of the next block |
712 | + elsif line =~ /^\[5D\]/ || line =~ /^\[9G\]/ | |
713 | + spherical = true | |
714 | + # Change the orbital type if necessary | |
715 | + gtos.each_with_index { | atom_gtos, atom_index| | |
716 | + atom_gtos.each { |gtoline| | |
717 | + if gtoline[0] >= 2 | |
718 | + gtoline[0] = -gtoline[0] # D5/F7/G9 | |
719 | + end | |
720 | + } | |
721 | + } | |
712 | 722 | elsif line =~ /^\[MO\]/ |
723 | + # Add shell info and primitive coefficients to molecule | |
724 | + gtos.each_with_index { | atom_gtos, atom_index| | |
725 | + atom_gtos.each { |gtoline| | |
726 | + sym = gtoline[0] | |
727 | + coeffs = gtoline[1] | |
728 | + nprimitives = coeffs.length / 2 | |
729 | + add_exp = gtoline[2] | |
730 | + ncomps += ncomp_hash[sym] | |
731 | + if !label | |
732 | + add_gaussian_orbital_shell(atom_index, sym, nprimitives, add_exp) | |
733 | + nprimitives.times { |prim| | |
734 | + add_gaussian_primitive_coefficients(coeffs[prim * 2], coeffs[prim * 2 + 1], 0.0) | |
735 | + } | |
736 | + end | |
737 | + } | |
738 | + } | |
713 | 739 | m = [] |
714 | 740 | idx_alpha = 1 # set_mo_coefficients() accepts 1-based index of MO |
715 | 741 | idx_beta = 1 |
@@ -767,8 +793,9 @@ class Molecule | ||
767 | 793 | end |
768 | 794 | break if i < ncomps # no MO info was found |
769 | 795 | end |
796 | + # TODO: reorder D, F, G coefficients for Molby order | |
770 | 797 | next |
771 | - end | |
798 | + end # end if | |
772 | 799 | end # end while |
773 | 800 | end # end catch |
774 | 801 | if errmsg |
@@ -780,27 +807,10 @@ class Molecule | ||
780 | 807 | |
781 | 808 | # Import the JANPA log and related molden files |
782 | 809 | # Files: inppath.{NAO.molden,CLPO.molden,janpa.log} |
783 | - # If mo_path is not nil, then clear existing mo info and load from mo_path | |
784 | - def sub_load_janpa_log(inppath, mo_path = nil) | |
810 | + # If inppath.spherical.molden is available, then clear existing mo info | |
811 | + # and load from it (i.e. use the basis set converted by molden2molden) | |
812 | + def sub_load_janpa_log(inppath) | |
785 | 813 | begin |
786 | - if mo_path | |
787 | - fp = File.open(mo_path, "rt") rescue fp = nil | |
788 | - if fp == nil | |
789 | - hide_progress_panel # Close if it is open | |
790 | - message_box("Cannot open MO molden file #{mo_path}: " + $!.to_s) | |
791 | - return false | |
792 | - end | |
793 | - @lineno = 0 | |
794 | - type = get_mo_info(:type) | |
795 | - alpha = get_mo_info(:alpha) | |
796 | - beta = get_mo_info(:beta) | |
797 | - clear_basis_set | |
798 | - set_mo_info(:type=>type, :alpha=>alpha, :beta=>beta) | |
799 | - # mol.@hf_type should be set before calling sub_load_molden | |
800 | - @hf_type = type | |
801 | - sub_load_molden(fp) | |
802 | - fp.close | |
803 | - end | |
804 | 814 | fp = File.open(inppath + ".janpa.log", "rt") rescue fp = nil |
805 | 815 | if fp == nil |
806 | 816 | hide_progress_panel # Close if it is open |
@@ -817,7 +827,30 @@ class Molecule | ||
817 | 827 | nao_infos = [] # index=atom_index, value=Hash with key "s", "px", "py" etc. |
818 | 828 | # nao_infos[index][key]: array of [nao_num, occupancy], in the reverse order of appearance |
819 | 829 | while line = getline.call |
820 | - if line =~ /^NAO \#/ | |
830 | + if line =~ /molden2molden: a conversion tool for MOLDEN/ | |
831 | + while line = getline.call | |
832 | + break if line =~ /^All done!/ | |
833 | + if line =~ /\.spherical\.molden/ | |
834 | + # The MOs are converted to spherical basis set | |
835 | + # Clear the existing MO and load *.spherical.molden | |
836 | + sname = inppath + ".spherical.molden" | |
837 | + fps = File.open(sname, "rt") rescue fps = nil | |
838 | + if fps != nil | |
839 | + print("Importing #{sname}.\n") | |
840 | + @lineno = 0 | |
841 | + type = get_mo_info(:type) | |
842 | + alpha = get_mo_info(:alpha) | |
843 | + beta = get_mo_info(:beta) | |
844 | + clear_basis_set | |
845 | + set_mo_info(:type=>type, :alpha=>alpha, :beta=>beta) | |
846 | + # mol.@hf_type should be set before calling sub_load_molden | |
847 | + @hf_type = type | |
848 | + sub_load_molden(fps) | |
849 | + fps.close | |
850 | + end | |
851 | + end | |
852 | + end | |
853 | + elsif line =~ /^NAO \#/ | |
821 | 854 | h["NAO"] = [] |
822 | 855 | while line = getline.call |
823 | 856 | break if line !~ /^\s*[1-9]/ |