Changeset 9fdbcc

Timestamp:

May 14, 2013, 12:46:47 PM (12 years ago)

Author:

Frederik Heber <heber@…>

Children:

28666f

Parents:

d13755

git-author:

Frederik Heber <heber@…> (03/27/13 11:10:25)

git-committer:

Frederik Heber <heber@…> (05/14/13 12:46:47)

Message:

Added use of switch in MPQCData to change between valence only and total charges.

File:

• src/bin/mpqc/mpqc.cc (modified) (6 diffs)

src/bin/mpqc/mpqc.cc

@@ -1180 +1180 @@
     delete[] ckptfile;
   }
+}
+
+static int getCoreElectrons(const int z)
+{
+  int n=0;
+  if (z > 2) n += 2;
+  if (z > 10) n += 8;
+  if (z > 18) n += 8;
+  if (z > 30) n += 10;
+  if (z > 36) n += 8;
+  if (z > 48) n += 10;
+  if (z > 54) n += 8;
+  return n;
 }
 
@@ -1640 +1653 @@
          data.charges.reserve(wfn->molecule()->natom());
          for (int iatom=0;iatom < wfn->molecule()->natom(); ++iatom) {
-           data.charges.push_back(wfn->molecule()->Z(iatom));
+           double charge = wfn->molecule()->Z(iatom);
+           if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly)
+             charge -= getCoreElectrons((int)charge);
+           data.charges.push_back(charge);
            std::vector<double> pos(3, 0.);
            for (int j=0;j<3;++j)
@@ -1714 +1730 @@
          SCVector3 r = min;
 
+         std::set<int> valence_indices;
+         if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly) {
+           // find valence orbitals
+           RefDiagSCMatrix evals = scf->eigenvalues();
+//           std::cout << "All Eigenvalues:" << std::endl;
+//           for(int i=0;i<wfn->oso_dimension(); ++i)
+//             std::cout << i << "th eigenvalue is " << evals(i) << std::endl;
+           std::set<double> evals_sorted;
+           for(int i=0;i<wfn->oso_dimension(); ++i)
+             evals_sorted.insert(evals(i));
+           std::set<double> evals_distances;
+           std::set<double>::const_iterator advancer = evals_sorted.begin();
+           std::set<double>::const_iterator iter = advancer++;
+           for(;advancer != evals_sorted.end(); ++advancer,++iter)
+             evals_distances.insert((*advancer)-(*iter));
+           const double largest_distance = *(evals_distances.rbegin());
+           ExEnv::out0()  << "Largest distance between EV is " << largest_distance << std::endl;
+           advancer = evals_sorted.begin();
+           iter = advancer++;
+           for(;advancer != evals_sorted.begin(); ++advancer,++iter)
+             if (fabs(fabs((*advancer)-(*iter)) - largest_distance) < 1e-10)
+               break;
+           assert( advancer != evals_sorted.begin() );
+           const double last_core_ev = (*iter);
+           ExEnv::out0()  << "Last core EV is " << last_core_ev << std::endl;
+           for(int i=0;i<wfn->oso_dimension(); ++i)
+             if (evals(i) > last_core_ev)
+               valence_indices.insert(i);
+//           {
+//             int i=0;
+//             std::cout << "Valence eigenvalues:" << std::endl;
+//             for (std::set<int>::const_iterator iter = valence_indices.begin();
+//                 iter != valence_indices.end(); ++iter)
+//               std::cout << i++ << "th eigenvalue is " << (*iter) << std::endl;
+//           }
+         } else {
+           // just insert all indices
+           for(int i=0;i<wfn->oso_dimension(); ++i)
+             valence_indices.insert(i);
+         }
+
+         // testing alternative routine from SCF::so_density()
+         RefSCMatrix oso_vector = scf->oso_eigenvectors();
+         RefSCMatrix vector = scf->so_to_orthog_so().t() * oso_vector;
+         oso_vector = 0;
+         RefSymmSCMatrix occ(scf->oso_dimension(), scf->basis_matrixkit());
+         occ.assign(0.0);
+         for (std::set<int>::const_iterator iter = valence_indices.begin();
+             iter != valence_indices.end(); ++iter) {
+           const int i = *iter;
+           occ(i,i) = scf->occupation(i);
+         }
+         RefSymmSCMatrix d2(scf->so_dimension(), scf->basis_matrixkit());
+         d2.assign(0.0);
+         d2.accumulate_transform(vector, occ);
+
          // taken from scf::density()
          RefSCMatrix nos
@@ -1735 +1807 @@
                    double tmp = 0.0;
                    for (int j=0; j<nbasis; ++j) {
-                       tmp += nos.get_element(j,i)*bs_values[j];
+                       tmp += d2(j,i)*bs_values[j]*bs_values[i];
                      }
-                   elec_density += nd.get_element(i)*tmp*tmp;
+                   elec_density += tmp;
                  }
                const double dens_at_r = elec_density * element_volume_conversion;
@@ -1764 +1836 @@
                integral_value += *diter;
            integral_value *= volume_element;
+           int n_electrons = scf->nelectron();
+           if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly)
+             n_electrons -= wfn->molecule()->n_core_electrons();
            const double normalization =
-               ((integral_value == 0) && (scf->nelectron() == 0)) ?
-                   1. : scf->nelectron()/integral_value;
+               ((integral_value == 0) && (n_electrons == 0)) ?
+                   1. : n_electrons/integral_value;
            std::cout << "Created " << data.sampled_grid.sampled_grid.size() << " grid points"
                << " with integral value of " << integral_value
-               << " against nelectron of " << scf->nelectron() << "." << std::endl;
+               << " against " << ((data.DoValenceOnly == MPQCData::DoSampleValenceOnly) ? "n_valence_electrons" : "n_electrons")
+               << " of " << n_electrons << "." << std::endl;
            for (std::vector<double>::iterator diter = data.sampled_grid.sampled_grid.begin();
                diter != data.sampled_grid.sampled_grid.end(); ++diter)
@@ -1907 +1983 @@
   MPQCData data(grid);
   data.DoLongrange = DoLongrange; // set whether we sample the density
+  data.DoValenceOnly = DoValenceOnly; // set whether we sample just valence electron and nuclei densities
 // now call work horse
   mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);