The electron density surface can be drawn by create_surface(:total_density).

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@555 a2be9bc6-48de-4e38-9406-05402d4bc13c

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index 373c808..11ca69f 100755 (executable)
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -11858,19 +11858,26 @@ MoleculeUpdateMCube(Molecule *mol, int idn)
                if (mc->radii == NULL)
                        return -2;  /*  Out of memory  */
                mc->nradii = mol->natoms;
-               memset(mc->radii, 0, sizeof(Double) * mc->nradii);
-               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)
-                               continue;  /*  This may happen when molecule is edited after setting up MO info  */
-                       mop = mobasep + sp->m_idx;
-                       for (iy = 0; iy < sp->ncomp; iy++) {
-                               dd = fabs(mop[iy]);
-                               if (dd > mc->radii[sp->a_idx])
-                                       mc->radii[sp->a_idx] = dd;
-                               if (dd > mopmax)
-                                       mopmax = dd;
+               if (mc->idn == mol->bset->nmos + 1) {
+                       /*  Total electron density  */
+                       for (ix = 0; ix < mol->natoms; ix++)
+                               mc->radii[ix] = 1.0;
+                       mopmax = 1.0;
+               } else {
+                       memset(mc->radii, 0, sizeof(Double) * mc->nradii);
+                       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)
+                                       continue;  /*  This may happen when molecule is edited after setting up MO info  */
+                               mop = mobasep + sp->m_idx;
+                               for (iy = 0; iy < sp->ncomp; iy++) {
+                                       dd = fabs(mop[iy]);
+                                       if (dd > mc->radii[sp->a_idx])
+                                               mc->radii[sp->a_idx] = dd;
+                                       if (dd > mopmax)
+                                               mopmax = dd;
+                               }
                        }
                }
                xmin = ymin = zmin = 1e10;
@@ -11950,7 +11957,33 @@ MoleculeUpdateMCube(Molecule *mol, int idn)
                                        n = (ix * ny + iy) * nz + iz;
                                        if (mc->dp[n] == DBL_MAX) {
                                                p.z = mc->origin.z + mc->dz * iz;
-                                               mc->dp[n] = sCalcMOPoint(mol, mol->bset, mc->idn, &p, tmp);
+                                               if (mc->idn == mol->bset->nmos + 1) {
+                                                       /*  Total electron density  */
+                                                       Int ne_alpha, ne_beta;
+                                                       mc->dp[n] = 0.0;
+                                                       ne_alpha = mol->bset->ne_alpha;
+                                                       ne_beta = mol->bset->ne_beta;
+                                                       if (mol->bset->rflag == 2 && ne_alpha < ne_beta) {
+                                                               /*  ROHF case: ensure ne_alpha >= ne_beta  */
+                                                               ne_beta = ne_alpha;
+                                                               ne_alpha = mol->bset->ne_beta;
+                                                       }
+                                                       for (sn = 1; sn <= ne_alpha; sn++) {
+                                                               dd = sCalcMOPoint(mol, mol->bset, sn, &p, tmp);
+                                                               dd = dd * dd;
+                                                               if (mol->bset->rflag != 0 && sn <= ne_beta)
+                                                                       dd *= 2;
+                                                               mc->dp[n] += dd;
+                                                       }
+                                                       if (mol->bset->rflag == 0) {
+                                                               for (sn = 1; sn <= ne_beta; sn++) {
+                                                                       dd = sCalcMOPoint(mol, mol->bset, sn + mol->bset->ncomps, &p, tmp);
+                                                                       mc->dp[n] += dd * dd;
+                                                               }
+                                                       }
+                                               } else {
+                                                       mc->dp[n] = sCalcMOPoint(mol, mol->bset, mc->idn, &p, tmp);
+                                               }
                                        }
                                }
                        }

diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c
index 197b383..13b66fc 100644 (file)
--- a/MolLib/Ruby_bind/ruby_bind.c
+++ b/MolLib/Ruby_bind/ruby_bind.c
@@ -10218,6 +10218,8 @@ s_Molecule_CreateSurface(int argc, VALUE *argv, VALUE self)
                rb_raise(rb_eMolbyError, "No MO information is given");
        if (nval == Qnil) {
                nmo = -1;
+       } else if (nval == ID2SYM(rb_intern("total_density"))) {
+               nmo = mol->bset->nmos + 1;
        } else {
                nmo = NUM2INT(rb_Integer(nval));
                if (nmo > mol->bset->nmos || nmo < -mol->bset->ncomps)

diff --git a/Scripts/commands.rb b/Scripts/commands.rb
index e36a60f..34bbd1f 100755 (executable)
--- a/Scripts/commands.rb
+++ b/Scripts/commands.rb
@@ -229,14 +229,12 @@ class Molecule
          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)
+         mult = (motype == "UHF" ? 2 : 1)
          mo_menu = ["1000 (-00.00000000)"]  #  Dummy entry; create later
          tabvals = []
          coeffs = nil
          a_idx_old = -1
+         occ = nil
          ncomps.times { |i|
            a_idx, s_idx, label = mol.get_gaussian_component_info(i)
                if a_idx_old != a_idx
@@ -297,15 +295,35 @@ class Molecule
          on_mo_action = lambda { |it|
            mo = it[:value]
                if mo_ao == 0
-                 if beta
+                 if motype == "UHF"
                    mo_index = (mo / 2) + (mo % 2 == 1 ? ncomps : 0) + 1
+                       if mo_index <= alpha || (mo_index > ncomps && mo_index <= ncomps + beta)
+                         occ_new = 1
+                       else
+                         occ_new = 0
+                       end
                  else
                    mo_index = mo + 1
+                       if mo_index <= alpha && mo_index <= beta
+                         occ_new = 1
+                       elsif mo_index <= alpha || mo_index <= beta
+                         occ_new = -1
+                       else
+                         occ_new = 0
+                       end
                  end
                else
                  mo_index = mo
+                 occ_new = occ
                end
                coeffs = nil
+               if occ_new != occ
+                 #  Set default color
+                 col = (occ_new == 0 ? 1 : (occ_new == -1 ? 2 : 0))
+                 set_value("color", col)
+                 h["color"] = col
+                 occ = occ_new
+               end
                item_with_tag("table")[:refresh] = true
            on_action.call(it)
          }
@@ -313,8 +331,8 @@ class Molecule
            if mo_ao != it[:value]
                  mo_ao = it[:value]
                  if mo_ao == 0
-                   mo_menu = (1..(ncomps * i)).map { |n|
-                 if beta
+                   mo_menu = (1..(ncomps * mult)).map { |n|
+                 if motype == "UHF"
                        i1 = (n - 1) / 2 + 1
                        i2 = n % 2
                        c1 = (i2 == 0 ? "B" : "A")
@@ -323,7 +341,13 @@ class Molecule
                        i1 = n
                        i2 = 1
                        c1 = ""
-                       c2 = (i1 > alpha ? "*" : "")
+                               if i1 > beta && i1 > alpha
+                                 c2 = "*"
+                               elsif i1 > beta || i1 > alpha
+                                 c2 = "S"
+                               else
+                                 c2 = ""
+                               end
                      end
                      en = mol.get_mo_energy(i1 + (i2 == 0 ? ncomps : 0))
                      sprintf("%d%s%s (%.8f)", i1, c1, c2, en)