OSDN Git Service

Show MO Surface dialog is implemented.
authortoshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Thu, 19 Jun 2014 05:50:27 +0000 (05:50 +0000)
committertoshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Thu, 19 Jun 2014 05:50:27 +0000 (05:50 +0000)
git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@551 a2be9bc6-48de-4e38-9406-05402d4bc13c

MolLib/Molecule.c
MolLib/Molecule.h
MolLib/Ruby_bind/ruby_bind.c
Scripts/commands.rb
Scripts/crystal.rb
Scripts/loadsave.rb
Scripts/molecule.rb
Scripts/uff.rb

index 222d747..74fc495 100755 (executable)
@@ -2444,7 +2444,85 @@ MoleculeAddGaussianPrimitiveCoefficients(Molecule *mol, Double exponent, Double
        return 0;
 }
 
-/*  Set MO coefficients for idx-th MO  */
+/*  Get the shell information from the component index  */
+/*  The outLabel must have space for at least 23 non-Null characters  */
+int
+MoleculeGetGaussianShellInfo(Molecule *mol, Int comp_idx, Int *outAtomIdx, char *outLabel, Int *outNprims)
+{
+       BasisSet *bset;
+       ShellInfo *shellp;
+       int si;
+       if (mol == NULL || (bset = mol->bset) == NULL)
+               return -1;  /*  No basis set info  */
+       if (comp_idx < 0 || comp_idx >= bset->ncomps)
+               return -2;  /*  Component index out of range  */
+       for (si = 0, shellp = bset->shells; si < bset->nshells; si++, shellp++) {
+               if (comp_idx >= shellp->ncomp) {
+                       comp_idx -= shellp->ncomp;
+                       continue;
+               } else {
+                       static const char *type_p = "xyz";
+                       static const char *type_d = "xxyyzzxyxzyz";
+                       static const char *type_d5[] = {"xy","yz","zz", "xz", "xx-yy"};
+                       static const char *type_f = "xxxyyyzzzxxyxxzxyyyyzxzzyzzxyz";
+                       static const char *type_f7[] = {"x3-3xy2", "x2z-y2z", "x(5z2-r2)", "z(5z2-3r2)", "y(5z2-r2)", "xyz", "3x2y-y3"};
+                       static const char *type_g[] = {"x4", "y4", "z4", "x3y", "x3z", "xy3", "y3z", "xz3", "yz3", "x2y2", "x2z2", "y2z2", "x2yz", "x2yz", "xyz2"};
+                       static const char *type_g9[] = {"x4+y4-6x2y2", "xz(x2-3y2)", "(x2-y2)(7z2-r2)", "xz(7z2-3r2)", "35z4-30z2r2+3r4", "yz(7z2-3r2)", "xy(7z2-r2)", "yz(3x2-y2)", "xy(x2-y2)"};
+                       *outAtomIdx = shellp->a_idx;
+                       *outNprims = shellp->nprim;
+                       switch (shellp->sym) {
+                               case kGTOType_S:
+                                       strcpy(outLabel, "S");
+                                       break;
+                               case kGTOType_P:
+                                       outLabel[0] = 'P';
+                                       outLabel[1] = type_p[comp_idx];
+                                       outLabel[2] = 0;
+                                       break;
+                               case kGTOType_SP:
+                                       if (comp_idx == 0)
+                                               strcpy(outLabel, "S");
+                                       else {
+                                               outLabel[0] = 'P';
+                                               outLabel[1] = type_p[comp_idx - 1];
+                                               outLabel[2] = 0;
+                                       }
+                                       break;
+                               case kGTOType_D:
+                                       outLabel[0] = 'D';
+                                       strncpy(outLabel + 1, type_d + comp_idx * 2, 2);
+                                       outLabel[3] = 0;
+                                       break;
+                               case kGTOType_D5:
+                                       outLabel[0] = 'D';
+                                       strcpy(outLabel + 1, type_d5[comp_idx]);
+                                       break;
+                               case kGTOType_F:
+                                       outLabel[0] = 'F';
+                                       strncpy(outLabel + 1, type_f + comp_idx * 3, 3);
+                                       outLabel[4] = 0;
+                               case kGTOType_F7:
+                                       outLabel[0] = 'F';
+                                       strcpy(outLabel + 1, type_f7[comp_idx]);
+                                       break;
+                               case kGTOType_G:
+                                       outLabel[0] = 'G';
+                                       strcpy(outLabel + 1, type_g[comp_idx]);
+                                       break;
+                               case kGTOType_G9:
+                                       outLabel[0] = 'G';
+                                       strcpy(outLabel + 1, type_g9[comp_idx]);
+                                       break;
+                               default:
+                                       return -3;  /*  Unsupported orbital type (internal error) */
+                       }
+                       return 0;
+               }
+       }
+       return -4;  /*  comp_idx out of range? (internal error)  */
+}
+
+/*  Set MO coefficients for idx-th MO (1-based)  */
 int
 MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Double *coeffs)
 {
@@ -2473,8 +2551,8 @@ MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Dou
                        bset->nmos = bset->ncomps;
                if (bset->nmos <= 0)
                        return -3;  /*  Bad or inconsistent number of MOs  */
-               bset->mo = (Double *)calloc(sizeof(Double), bset->nmos * bset->ncomps);
-               bset->moenergies = (Double *)calloc(sizeof(Double), bset->nmos);
+               bset->mo = (Double *)calloc(sizeof(Double), (bset->nmos + 1) * bset->ncomps);
+               bset->moenergies = (Double *)calloc(sizeof(Double), bset->nmos + 1);
                if (bset->mo == NULL || bset->moenergies == NULL) {
                        if (bset->mo != NULL)
                                free(bset->mo);
@@ -2486,8 +2564,14 @@ MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Dou
                        return -2;  /*  Low memory  */
                }
        }
-       if (idx < 0 || idx >= bset->nmos)
+       if (idx < 0)
+               idx = -idx + bset->ncomps;
+       if (idx < 0 || idx > bset->nmos)
                return -4;  /*  Bad MO index  */
+       if (idx == 0)
+               idx = bset->nmos;  /*  Arbitrary vector  */
+       else
+               idx--;
        if (energy != -1000000)
                bset->moenergies[idx] = energy;
        if (ncomps < bset->ncomps)
@@ -2502,7 +2586,7 @@ MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Dou
        return 0;
 }
 
-/*  Get MO coefficients for idx-th MO  */
+/*  Get MO coefficients for idx-th MO (1-based)  */
 /*  Caution: *ncoeffs and *coeffs should be valid _before_ calling this function, i.e.  */
 /*  *ncoeffs = 0 && *coeffs = NULL or *coeffs is a valid memory pointer and *ncoeffs  */
 /*  properly designates the memory size as an array of Doubles.  */
@@ -2515,8 +2599,14 @@ MoleculeGetMOCoefficients(Molecule *mol, Int idx, Double *energy, Int *ncoeffs,
        bset = mol->bset;
        if (bset == NULL || bset->ncomps <= 0)
                return -2;  /*  No basis set info  */
-       if (idx < 0 || idx >= bset->nmos)
+       if (idx < 0)
+               idx = -idx + bset->ncomps;
+       if (idx < 0 || idx > bset->nmos)
                return -3;  /*  MO index out of range  */
+       if (idx == 0)
+               idx = bset->nmos;  /*  Arbitrary vector  */
+       else
+               idx--;
        if (energy != NULL)
                *energy = bset->moenergies[idx];
        if (ncoeffs != NULL && coeffs != NULL) {
@@ -3277,7 +3367,7 @@ MoleculeLoadGamessDatFile(Molecule *mol, const char *fname, char **errbuf)
                                }
                                if (k < mol->bset->ncomps)
                                        continue;
-                               j = MoleculeSetMOCoefficients(mol, i, -1000000, k, coeffs);
+                               j = MoleculeSetMOCoefficients(mol, i + 1, -1000000, k, coeffs);
                                if (j != 0) {
                                        s_append_asprintf(errbuf, "Line %d: cannot set coefficients for MO %d", lineNumber, i + 1);
                                        free(coeffs);
@@ -10787,7 +10877,7 @@ MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Doub
 #pragma mark ====== MO calculation ======
 
 /*  Calculate an MO value for a single point.  */
-/*  Index is the MO number (1-based)  */
+/*  Index is the MO number (1-based); 0 denotes "arbitrary vector"  */
 /*  tmp is an array of (natoms * 4) atoms, and used to store dr and |dr|^2 for each atom.  */
 static Double
 sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
@@ -10797,6 +10887,8 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do
        Double val, tval, *cnp, *tmpp, *mobasep, *mop;
        Int i, j;
        /*  Cache dr and |dr|^2  */
+       if (index == 0)
+               index = bset->nmos + 1;
        for (i = 0; i < mp->natoms; i++) {
                Vector r;
                r = ATOM_AT_INDEX(mp->atoms, i)->r;
@@ -10908,7 +11000,7 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do
 }
 
 /*  Calculate one MO. The input vectors are angstrom unit (changed from bohr unit: 20140520)  */
-/*  mono is the MO number (1-based)  */
+/*  mono is the MO number (1-based); 0 denotes "arbitrary vector" */
 int
 MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, const Vector *dyp, const Vector *dzp, Int nx, Int ny, Int nz, int (*callback)(double progress, void *ref), void *ref)
 {
@@ -11476,7 +11568,7 @@ MoleculeUpdateMCube(Molecule *mol, int idn)
                return -1;  /*  Number of atoms is smaller than expected  */
 
        mc = mol->mcube;
-       if (idn > 0) {
+       if (idn >= 0) {
                ShellInfo *sp;
                Double *mobasep, *mop, mopmax;
                Double xmin, xmax, ymin, ymax, zmin, zmax;
@@ -11493,7 +11585,7 @@ MoleculeUpdateMCube(Molecule *mol, int idn)
                        return -2;  /*  Out of memory  */
                mc->nradii = mol->natoms;
                memset(mc->radii, 0, sizeof(Double) * mc->nradii);
-               mobasep = mol->bset->mo + (mc->idn - 1) * mol->bset->ncomps;
+               mobasep = mol->bset->mo + (mc->idn == 0 ? mol->bset->nmos : mc->idn - 1) * mol->bset->ncomps;
                mopmax = 0.0;
                for (ix = 0, sp = mol->bset->shells; ix < mol->bset->nshells; ix++, sp++) {
                        if (sp->a_idx >= mol->natoms)
@@ -11782,8 +11874,7 @@ MoleculeUpdateMCube(Molecule *mol, int idn)
                        /*  Less than 3 points: no triangles  */
                        if (mc->c[sn].ntriangles > 0)
                                mc->c[sn].triangles[0] = -1;  /*  End mark  */
-                       retval = 0;
-                       goto end;
+                       continue;
                }
                
                /*  Create triangle table  */
index 18cac8f..5fb9d0b 100755 (executable)
@@ -239,6 +239,7 @@ typedef struct BasisSet {
        Int ncomps;          /*  Number of AO components; equal to sum of shells[i].ncomp  */
        Int nmos;            /*  Number of MOs; equal to ncomps if close shell, ncomps*2 if open shell */
        Double *mo;          /*  MO matrix (mo[i][j] represents the j-th AO coefficient for the i-th MO)  */
+                                                /*  Memory are allocated for (2*nmos+1) entries; the last entry is for displaying arbitrary vector  */
        Double *moenergies;  /*  MO energies  */
        Double *scfdensities; /*  SCF densities; lower triangle of a symmetric matrix (size nmos*(nmos+1)/2)  */
        Int ncubes;          /*  Number of calculated MOs  */
@@ -409,6 +410,7 @@ void MoleculeExchange(Molecule *mp1, Molecule *mp2);
 
 int MoleculeAddGaussianOrbitalShell(Molecule *mol, Int sym, Int nprims, Int a_idx);
 int MoleculeAddGaussianPrimitiveCoefficients(Molecule *mol, Double exponent, Double contraction, Double contraction_sp);
+int MoleculeGetGaussianShellInfo(Molecule *mol, Int comp_idx, Int *outAtomIdx, char *outLabel, Int *outNprims);
 int MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Double *coeffs);
 int MoleculeGetMOCoefficients(Molecule *mol, Int idx, Double *energy, Int *ncoeffs, Double **coeffs);
 int MoleculeAllocateBasisSetRecord(Molecule *mol, Int rflag, Int ne_alpha, Int ne_beta);
index 60ca518..5762a17 100644 (file)
@@ -10132,8 +10132,9 @@ s_Molecule_Cubegen(int argc, VALUE *argv, VALUE self)
  *  call-seq:
  *     create_surface(mo, attr = nil)
  *
- *  Create a MO surface. The argument mo is the MO index (1-based); if mo is -1,
- *  then the attributes of the current surface are modified.
+ *  Create a MO surface. The argument mo is the MO index (1-based); if mo is negative,
+ *  then it denotes the beta orbital.
+ *  If mo is nil, then the attributes of the current surface are modified.
  *  Attributes:
  *    :npoints : the approximate number of grid points
  *    :expand  : the scale factor to expand/shrink the display box size for each atom,
@@ -10154,11 +10155,17 @@ s_Molecule_CreateSurface(int argc, VALUE *argv, VALUE self)
        Double d[4];
     Data_Get_Struct(self, Molecule, mol);
        rb_scan_args(argc, argv, "11", &nval, &hval);
-       nmo = NUM2INT(rb_Integer(nval));
        if (mol->bset == NULL)
                rb_raise(rb_eMolbyError, "No MO information is given");
-       if ((nmo <= 0 && nmo != -1) || nmo > mol->bset->nmos)
-               rb_raise(rb_eMolbyError, "MO index (%d) is out of range; should be 1..%d", nmo, mol->bset->nmos);
+       if (nval == Qnil) {
+               nmo = -1;
+       } else {
+               nmo = NUM2INT(rb_Integer(nval));
+               if (nmo > mol->bset->nmos || nmo < -mol->bset->ncomps)
+                       rb_raise(rb_eMolbyError, "MO index (%d) is out of range; should be 1..%d (or -1..-%d for beta orbitals); (0 is acceptable as arbitrary vector)", nmo, mol->bset->nmos, mol->bset->ncomps);
+               if (nmo < 0)
+                       nmo = -nmo + mol->bset->ncomps;
+       }
        if (hval != Qnil && (aval = rb_hash_aref(hval, ID2SYM(rb_intern("npoints")))) != Qnil) {
                npoints = NUM2INT(rb_Integer(aval));
                need_recalc = 1;
@@ -10219,7 +10226,7 @@ static VALUE
 s_Molecule_SetSurfaceAttr(VALUE self, VALUE hval)
 {
        VALUE args[2];
-       args[0] = INT2FIX(-1);
+       args[0] = Qnil;
        args[1] = hval;
        return s_Molecule_CreateSurface(2, args, self);
 }
@@ -10340,6 +10347,30 @@ s_Molecule_AddGaussianPrimitiveCoefficients(VALUE self, VALUE expval, VALUE cval
 
 /*
  *  call-seq:
+ *     get_gaussian_shell_info(comp_index) -> [atom_index, orbital_description, no_of_primitives]
+ *
+ *  Get the Gaussian shell information for the given MO coefficient index.
+ */
+static VALUE
+s_Molecule_GetGaussianShellInfo(VALUE self, VALUE cval)
+{
+       Molecule *mol;
+       Int n, c, atom_idx, nprims;
+       char label[32];
+    Data_Get_Struct(self, Molecule, mol);
+       if (mol->bset == NULL)
+               return Qnil;
+       c = NUM2INT(rb_Integer(cval));
+       if (c < 0 || c >= mol->bset->ncomps)
+               rb_raise(rb_eMolbyError, "The component index (%d) is out of range (should be 0..%d)", c, mol->bset->ncomps - 1);
+       n = MoleculeGetGaussianShellInfo(mol, c, &atom_idx, label, &nprims);
+       if (n != 0)
+               rb_raise(rb_eMolbyError, "Cannot get the shell info for component index (%d)", c);
+       return rb_ary_new3(3, INT2NUM(atom_idx), rb_str_new2(label), INT2NUM(nprims));
+}
+
+/*
+ *  call-seq:
  *     clear_mo_coefficients
  *
  *  Clear the existing MO coefficients.
@@ -10367,8 +10398,9 @@ s_Molecule_ClearMOCoefficients(VALUE self)
  *  call-seq:
  *     set_mo_coefficients(idx, energy, coefficients)
  *
- *  To be used internally. Add a MO coefficients. Idx is the MO index (for open shell system, 
- *  beta MOs comes after all alpha MOs), energy is the MO energy, coefficients is an array
+ *  To be used internally. Add a MO coefficients. Idx is the MO index (1-based; for open shell system, 
+ *  beta MOs comes after all alpha MOs; alternatively, the beta MO can be specified as a negative MO number)
+ *  Energy is the MO energy, and coefficients is an array
  *  of MO coefficients.
  */
 static VALUE
@@ -10390,7 +10422,7 @@ s_Molecule_SetMOCoefficients(VALUE self, VALUE ival, VALUE eval, VALUE aval)
        }
        for (i = 0; i < ncomps; i++)
                coeffs[i] = NUM2DBL(rb_Float(RARRAY_PTR(aval)[i]));
-       i = MoleculeSetMOCoefficients(mol, idx, energy, ncomps, coeffs);
+       i = MoleculeSetMOCoefficients(mol, idx, energy, ncomps, coeffs); /* Negative (beta orbitals) or zero (arbitrary vector) idx is allowed */
 end:
        if (i == -1)
                rb_raise(rb_eMolbyError, "Molecule is emptry");
@@ -10399,7 +10431,7 @@ end:
        else if (i == -3)
                rb_raise(rb_eMolbyError, "Bad or inconsistent number of MOs");
        else if (i == -4)
-               rb_raise(rb_eMolbyError, "Bad MO index");
+               rb_raise(rb_eMolbyError, "Bad MO index (%d)", idx);
        else if (i == -5)
                rb_raise(rb_eMolbyError, "Insufficient number of coefficients are given");
        else if (i != 0)
@@ -10411,7 +10443,7 @@ end:
  *  call-seq:
  *     get_mo_coefficients(idx)
  *
- *  To be used internally. Get an array of MO coefficients for the given MO index (0-based).
+ *  To be used internally. Get an array of MO coefficients for the given MO index (1-based).
  */
 static VALUE
 s_Molecule_GetMOCoefficients(VALUE self, VALUE ival)
@@ -10443,7 +10475,7 @@ s_Molecule_GetMOCoefficients(VALUE self, VALUE ival)
  *  call-seq:
  *     get_mo_energy(idx)
  *
- *  To be used internally. Get the MO energy for the given MO index (0-based).
+ *  To be used internally. Get the MO energy for the given MO index (1-based).
  */
 static VALUE
 s_Molecule_GetMOEnergy(VALUE self, VALUE ival)
@@ -10463,7 +10495,7 @@ s_Molecule_GetMOEnergy(VALUE self, VALUE ival)
        return rb_float_new(energy);
 }
 
-static VALUE sTypeSym, sAlphaSym, sBetaSym;
+static VALUE sTypeSym, sAlphaSym, sBetaSym, sNcompsSym;
 
 static inline void
 s_InitMOInfoKeys(void)
@@ -10472,6 +10504,7 @@ s_InitMOInfoKeys(void)
                sTypeSym = ID2SYM(rb_intern("type"));
                sAlphaSym = ID2SYM(rb_intern("alpha"));
                sBetaSym = ID2SYM(rb_intern("beta"));
+               sNcompsSym = ID2SYM(rb_intern("ncomps"));
        }
 }
 
@@ -10529,6 +10562,7 @@ s_Molecule_SetMOInfo(VALUE self, VALUE hval)
  *     get_mo_info(key)
  *
  *  Get the MO info. The key is as described in set_mo_info.
+ *  Read-only: :ncomps the number of components (and the number of MOs)
  */
 static VALUE
 s_Molecule_GetMOInfo(VALUE self, VALUE kval)
@@ -10536,7 +10570,7 @@ s_Molecule_GetMOInfo(VALUE self, VALUE kval)
        Molecule *mol;
     Data_Get_Struct(self, Molecule, mol);
        if (mol->bset == NULL)
-               rb_raise(rb_eMolbyError, "No MO information is defined");
+               return Qnil;
        if (kval == sTypeSym) {
                switch (mol->bset->rflag) {
                        case 0: return rb_str_new2("UHF");
@@ -10548,14 +10582,20 @@ s_Molecule_GetMOInfo(VALUE self, VALUE kval)
                return INT2NUM(mol->bset->ne_alpha);
        } else if (kval == sBetaSym) {
                return INT2NUM(mol->bset->ne_beta);
-       } else return Qnil;
+       } else if (kval == sNcompsSym) {
+               return INT2NUM(mol->bset->ncomps);
+       } else {
+               kval = rb_inspect(kval);
+               rb_raise(rb_eMolbyError, "Unknown MO info key: %s", StringValuePtr(kval));
+               return Qnil;  /*  Does not reach here  */
+       }
 }
 
 /*
  *  call-seq:
  *     mo_type
  *
- *  Returns either "RHF", "UHF", or "ROHF". If no MO info is present, raises exception.
+ *  Returns either "RHF", "UHF", or "ROHF". If no MO info is present, returns nil.
  */
 static VALUE
 s_Molecule_MOType(VALUE self)
@@ -11388,6 +11428,7 @@ Init_Molby(void)
        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_primitive_coefficients", s_Molecule_AddGaussianPrimitiveCoefficients, 3);
+       rb_define_method(rb_cMolecule, "get_gaussian_shell_info", s_Molecule_GetGaussianShellInfo, 1);
        rb_define_method(rb_cMolecule, "mo_type", s_Molecule_MOType, 0);
        rb_define_method(rb_cMolecule, "clear_mo_coefficients", s_Molecule_ClearMOCoefficients, 0);
        rb_define_method(rb_cMolecule, "set_mo_coefficients", s_Molecule_SetMOCoefficients, 3);
index 8ef7432..db841f5 100755 (executable)
@@ -130,9 +130,9 @@ class Molecule
            sprintf("%.8f", mol.get_property(col - 1, row)) rescue ""
          end
     }
-    Dialog.new("Extra Props:" + mol.name, nil, nil, :resizable=>true, :has_close_box=>true) {
+    mol.open_auxiliary_window("Extra Props", nil, nil, :resizable=>true, :has_close_box=>true) {
          names = nil
-         on_document_modified = lambda { |m|
+         @on_document_modified = lambda { |m|
            if (n = m.property_names) != names
                  names = n
                  col = [["frame", 40]] + names.map { |nn| [nn, 120] }
@@ -147,8 +147,8 @@ class Molecule
            :flex=>[0,0,0,0,1,1]
          )
          set_min_size(320, 200)
-         listen(mol, "documentModified", on_document_modified)
-         listen(mol, "documentWillClose", lambda { |m| close } )
+         listen(mol, "documentModified", on_document_modified)
+         listen(mol, "documentWillClose", lambda { |m| close } )
          on_document_modified.call(mol)
          show
        }
@@ -204,7 +204,7 @@ class Molecule
          }
          layout(1,
                item(:text, :title=>"Energy =                 ", :tag=>"energy"),
-               item(:view, :frame=>[0, 0, 100, 80], :tag=>"graph", :on_paint=>draw_graph, :flex=>[0,0,0,0,1,1]),
+               item(:view, :frame=>[0, 0, 320, 240], :tag=>"graph", :on_paint=>draw_graph, :flex=>[0,0,0,0,1,1]),
        #       item(:button, :title=>"Update", :action=>doc_modified, :align=>:center, :flex=>[1,1,1,0,0,0]),
        #       item(:button, :title=>"Close", :action=>proc { hide }, :align=>:center, :flex=>[1,1,1,0,0,0]),
                :flex=>[0,0,0,0,1,1]
@@ -219,6 +219,167 @@ class Molecule
        self
   end
   
+  def cmd_create_surface
+    mol = self
+       mol.open_auxiliary_window("MO Surface", nil, nil, :resizable=>true, :has_close_box=>true) {
+         motype = mol.get_mo_info(:type)
+         alpha = mol.get_mo_info(:alpha)
+         beta = mol.get_mo_info(:beta)
+         ncomps = mol.get_mo_info(:ncomps)
+         mo_index = 1
+         mo_ao = nil
+         coltable = [[0,0,1], [1,0,0], [0,1,0], [1,1,0], [0,1,1], [1,0,1], [0,0,0]]
+         if (motype == "RHF")
+           beta = nil
+         end
+         i = (beta ? 2 : 1)
+         mo_menu = []  #  Create later
+         tabvals = []
+         coeffs = nil
+         a_idx_old = -1
+         ncomps.times { |i|
+           a_idx, label, nprims = mol.get_gaussian_shell_info(i)
+               if a_idx_old != a_idx
+                 a_idx_old = a_idx
+                 a = a_idx.to_s
+                 n = mol.atoms[a_idx].name
+               else
+                 a = n = ""
+               end
+               tabvals.push([a, n, label, a_idx])
+         }
+         on_get_value = lambda { |it, row, col|
+           if col < 3
+                 tabvals[row][col]
+               else
+                 if coeffs == nil
+                   if mo_ao == 0
+                     coeffs = mol.get_mo_coefficients(mo_index)
+                   else
+                     coeffs = (0...ncomps).map { |i| (i == mo_index ? 1.0 : 0.0) }
+                       end
+                 end
+                 sprintf("%.6f", coeffs[row])
+               end
+         }
+         h = {"mo"=>nil, "color"=>nil, "opacity"=>nil, "threshold"=>nil, "expand"=>nil, "grid"=>nil}
+         should_update = true
+         on_action = lambda { |it|
+           should_update = false
+           h.each_key { |key|
+                 val = value(key)
+                 if val && h[key] != val
+                   should_update = true
+                   break
+                 end
+               }
+               item_with_tag("update")[:enabled] = should_update
+         }
+         on_mo_action = lambda { |it|
+           mo = it[:value]
+               if mo_ao == 0
+                 if beta
+                   mo_index = (mo / 2) + (mo % 2 == 1 ? ncomps : 0) + 1
+                 else
+                   mo_index = mo + 1
+                 end
+               else
+                 mo_index = mo
+               end
+               coeffs = nil
+               item_with_tag("table")[:refresh] = true
+           on_action.call(it)
+         }
+         on_set_action = lambda { |it|
+           if mo_ao != it[:value]
+                 mo_ao = it[:value]
+                 if mo_ao == 0
+                   mo_menu = (1..(ncomps * i)).map { |n|
+                 if beta
+                       i1 = (n - 1) / 2 + 1
+                       i2 = n % 2
+                       c1 = (i2 == 0 ? "B" : "A")
+                       c2 = (i1 > (i2 == 0 ? beta : alpha) ? "*" : "")
+                     else
+                       i1 = n
+                       i2 = 1
+                       c1 = ""
+                       c2 = (i1 >= alpha ? "*" : "")
+                     end
+                     en = mol.get_mo_energy(i1 + (i2 == 0 ? ncomps : 0))
+                     sprintf("%d%s%s (%.8f)", i1, c1, c2, en)
+                       }
+                 else
+                   mo_menu = []
+                   ncomps.times { |i|
+                         mo_menu[i] = sprintf("AO%d: %s (%s)", i + 1, tabvals[i][2], mol.atoms[tabvals[i][3]].name)
+                       }
+                 end
+                 it0 = item_with_tag("mo")
+                 it0[:subitems] = mo_menu
+                 it0[:value] = 0
+                 on_mo_action.call(it0)
+               end
+         }
+         on_update = lambda { |it|
+           h.each_key { |key|
+                 h[key] = value(key)
+               }
+               opac = h["opacity"].to_f
+               color = coltable[h["color"]] + [opac]
+               thres = h["threshold"].to_f
+               thres = 0.001 if thres >= 0.0 && thres < 0.001
+               thres = -0.001 if thres <= 0.0 && thres > -0.001
+               expand = h["expand"].to_f
+               expand = 0.01 if expand < 0.01
+               expand = 10.0 if expand > 10.0
+               grid = h["grid"].to_i
+               if grid > 10000000
+                 grid = 10000000
+               end
+               if mo_ao == 0
+                 idx = mo_index
+               else
+                 idx = 0
+                 mol.set_mo_coefficients(0, 0.0, coeffs)
+               end
+               mol.create_surface(idx, :npoints=>grid, :color=>color, :thres=>thres, :expand=>expand)
+               on_action.call(it)
+         }
+         layout(1,
+           layout(2,
+                 item(:text, :title=>"Orbital Set"),
+                 item(:popup, :tag=>"mo_ao", :subitems=>["Molecular Orbitals", "Atomic Orbitals"], :action=>on_set_action),
+             item(:text, :title=>"Select"),
+             item(:popup, :tag=>"mo", :subitems=>mo_menu, :action=>on_mo_action)),
+               layout(4,
+                 item(:text, :title=>"Color"),
+                 item(:popup, :tag=>"color", :subitems=>["blue", "red", "green", "yellow", "cyan", "magenta", "black"], :action=>on_action),
+                 item(:text, :title=>"Opacity"),
+                 item(:textfield, :tag=>"opacity", :width=>80, :value=>"0.6", :action=>on_action),
+                 item(:text, :title=>"Threshold"),
+                 item(:textfield, :tag=>"threshold", :width=>80, :value=>"0.05", :action=>on_action),
+                 item(:text, :title=>"Box Limit"),
+                 item(:textfield, :tag=>"expand", :width=>80, :value=>"1.0", :action=>on_action)),
+               layout(2,
+                 item(:text, :title=>"Number of Grid Points"),
+                 item(:textfield, :tag=>"grid", :width=>120, :value=>"512000", :action=>on_action)),
+           item(:table, :width=>300, :height=>300, :tag=>"table",
+                 :columns=>[["Atom", 60], ["Name", 60], ["Label", 60], ["Coeff", 120]],
+                 :on_count=> lambda { |it| tabvals.count },
+                 :on_get_value=>on_get_value,
+                 :flex=>[0,0,0,0,1,1]),
+               item(:button, :tag=>"update", :title=>"Update", :action=>on_update, :flex=>[0,1,0,0,0,0]),
+               :flex=>[0,0,0,0,1,1]
+         )
+         on_set_action.call(item_with_tag("mo_ao"))
+         size = self.size
+         set_min_size(size[0], 250)
+         item_with_tag("table")[:refresh] = true
+         show
+    }
+  end
+  
   #  DEBUG
   def cmd_test
     $test_dialog = Dialog.new("Test") { item(:text, :title=>"test"); show }
@@ -255,5 +416,6 @@ register_menu("Delete Frames...", :cmd_delete_frames, lambda { |m| m && m.nframe
 register_menu("Reverse Frames...", :cmd_reverse_frames, lambda { |m| m && m.nframes > 1 } )
 register_menu("", "")
 register_menu("Show Energy Window...", :cmd_show_energy, lambda { |m| m && m.property_names.include?("energy") } )
+register_menu("Show MO Surface...", :cmd_create_surface, lambda { |m| m && m.get_mo_info(:type) != nil } )
 #register_menu("cmd test", :cmd_test)
 
index 12c92a3..613c7c5 100755 (executable)
@@ -669,7 +669,7 @@ end
 def cmd_plane
   plane_settings = @plane_settings || Hash.new
   mol = self
-  h = Dialog.new("Best-Fit Planes: " + mol.name, "Close", nil) {
+  mol.open_auxiliary_window("Best-Fit Planes", "Close", nil) {
     refresh_proc = lambda { |it|
       n = it[:tag][/\d/].to_i
       g = plane_settings["group#{n}"]
@@ -1034,7 +1034,7 @@ end
 
 def cmd_bond_angle_with_sigma
   mol = self
-  Dialog.new("Bond & Angle with Sigma:" + self.name, nil, nil) {
+  mol.open_auxiliary_window("Bond & Angle with Sigma", nil, nil) {
     values = []
        clicked = []
        sel = mol.selection
@@ -1127,7 +1127,7 @@ def cmd_bond_angle_with_sigma
                  :enabled=>false),
                [item(:button, :title=>"Close", :action=>lambda { |item| hide }), {:align=>:right}])
        )
-       on_document_modified = lambda { |*d|
+       @on_document_modified = lambda { |*d|
          newsel = mol.selection - sel
          sel = mol.selection
          row = (item_with_tag("table")[:selection] || [nil])[-1]
@@ -1159,9 +1159,9 @@ def cmd_bond_angle_with_sigma
                item_with_tag("table")[:refresh] = true
          end
        }
-    listen(mol, "documentModified", on_document_modified)
-       listen(mol, "documentWillClose", lambda { |*d| hide } )
-       on_document_modified.call
+#   listen(mol, "documentModified", on_document_modified)
+#      listen(mol, "documentWillClose", lambda { |*d| hide } )
+       @on_document_modified.call
        show
   }
 end
@@ -1441,7 +1441,7 @@ def cmd_show_ortep
          ["H", "all", "Open", "black"]
        ]
   }
-  Dialog.new("Show ORTEP:" + mol.name, nil, nil, :resizable=>true, :has_close_box=>true) {
+  mol.open_auxiliary_window("Show ORTEP", nil, nil, :resizable=>true, :has_close_box=>true) {
     tepview = nil  #  Forward declaration
        tab = "Atoms"
        columns = {
@@ -1686,7 +1686,7 @@ def cmd_show_ortep
            filecopy("#{tmp}/TEP.IN", fname)
       end
        }
-       on_document_modified = lambda { |m|
+       @on_document_modified = lambda { |m|
        }
        #  Close handler (called when the close box is pressed or the document is closed)
        @on_close = lambda { |*d|
@@ -1781,9 +1781,9 @@ def cmd_show_ortep
        )
        set_min_size(480, 200)
        tepview = item_with_tag("ortep")
-       listen(mol, "documentModified", on_document_modified)
-       listen(mol, "documentWillClose", lambda { |m| close } )
-       on_document_modified.call(nil)
+#      listen(mol, "documentModified", on_document_modified)
+#      listen(mol, "documentWillClose", lambda { |m| close } )
+       @on_document_modified.call(nil)
     on_update_ortep.call(nil)
        show
   }
index db3eb00..c4dafcf 100755 (executable)
@@ -220,7 +220,7 @@ class Molecule
                        }
                end
                mo.each_with_index { |m, i|
-                       idx = Integer(mo_labels[i]) - 1
+                       idx = Integer(mo_labels[i])
                        set_mo_coefficients(idx + (alpha ? 0 : ncomps), Float(mo_energies[i]), m)
                #       if mo_labels[i] % 8 == 1
                #               puts "set_mo_coefficients #{idx}, #{mo_energies[i]}, [#{m[0]}, ..., #{m[-1]}]"
index 7116bca..d536ccd 100755 (executable)
@@ -757,7 +757,6 @@ class Molecule
   end
   
   def close_active_auxiliary_window
-    puts @aux_windows.inspect
     if @aux_windows
          @aux_windows.each_value { |d|
            puts "#{d}: #{d.active?}"
index e0477fa..9f11b8d 100644 (file)
@@ -279,7 +279,7 @@ def guess_uff_parameters
   arena.prepare(true)
   xatoms = IntGroup[]
   xbonds = xangles = xdihedrals = xfragments = []
-  h = Dialog.new("Guess UFF Parameters: #{mol.name}", nil, nil, :resizable=>true) {
+  mol.open_auxiliary_window("Guess UFF Parameters", nil, nil, :resizable=>true) {
     update_xatoms = lambda {
       xatoms = mol.atom_group { |ap| !exclude.member?(ap.atomic_number) }
       xfragments = mol.fragments(xatoms)
@@ -1091,8 +1091,9 @@ def guess_uff_parameters
     size = self.size
     set_min_size(size)
     set_size(size[0] + 100, size[1] + 50);
-    listen(mol, "documentModified", lambda { |d| update_xatoms.call; update_selection.call })
-    listen(mol, "documentWillClose", lambda { |d| hide })
+       @on_document_modified = lambda { |d| update_xatoms.call; update_selection.call }
+    # listen(mol, "documentModified", lambda { |d| update_xatoms.call; update_selection.call })
+    # listen(mol, "documentWillClose", lambda { |d| hide })
     update_xatoms.call
     guess_uff_types.call(mol.all)
     show