Changeset 62dabe

Timestamp:

Jul 6, 2012, 5:43:58 PM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

aa42a9

Parents:

d11832

Message:

Extracted Option parsing and setting into extra functions.

File:

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

src/bin/mpqc/mpqc.cc

@@ -49 +49 @@
 #include <fstream>
 
+#include <boost/bind.hpp>
+#include <boost/function.hpp>
+
 #include <scconfig.h>
 #ifdef HAVE_SSTREAM
@@ -178 +181 @@
 }
 
-int
-try_main(int argc, char *argv[])
+double EvaluateDensity(
+    SCVector3 &r,
+    Ref<Integral> &intgrl,
+    GaussianBasisSet::ValueData &vdat,
+    Ref<Wavefunction> &wfn);
+
+/** Places all known options into \a options and parses them from argc,argv.
+ *
+ * \param options options structure
+ * \param argc argument count
+ * \param argv argument array
+ * \return return value by GetLongOpt::parse() function
+ */
+int ParseOptions(GetLongOpt &options, int argc, char **argv)
 {
-  //trash_stack();
-
-  KeyValValueboolean truevalue(1), falsevalue(0);
-  int i;
-  const char *devnull = "/dev/null";
-  atexit(clean_up);
-
-#ifdef HAVE_FEENABLEEXCEPT
-  // this uses a glibc extension to trap on individual exceptions
-# ifdef FE_DIVBYZERO
-  feenableexcept(FE_DIVBYZERO);
-# endif
-# ifdef FE_INVALID
-  feenableexcept(FE_INVALID);
-# endif
-# ifdef FE_OVERFLOW
-  feenableexcept(FE_OVERFLOW);
-# endif
-#endif
-
-#ifdef HAVE_FEDISABLEEXCEPT
-  // this uses a glibc extension to not trap on individual exceptions
-# ifdef FE_UNDERFLOW
-  fedisableexcept(FE_UNDERFLOW);
-# endif
-# ifdef FE_INEXACT
-  fedisableexcept(FE_INEXACT);
-# endif
-#endif
-
-#if defined(HAVE_SETRLIMIT)
-  struct rlimit rlim;
-  rlim.rlim_cur = 0;
-  rlim.rlim_max = 0;
-  setrlimit(RLIMIT_CORE,&rlim);
-#endif
-
-  ExEnv::init(argc, argv);
-
-  Ref<MessageGrp> grp;
-#if defined(HAVE_MPI) && defined(ALWAYS_USE_MPI)
-  grp = new MPIMessageGrp(&argc, &argv);
-#endif
-
-  // parse commandline options
-  GetLongOpt options;
-
   options.usage("[options] [filename]");
   options.enroll("f", GetLongOpt::MandatoryValue,
@@ -254 +222 @@
   options.enroll("d", GetLongOpt::NoValue, "debug", 0);
   options.enroll("h", GetLongOpt::NoValue, "print this message", 0);
-  options.enroll("cca-path", GetLongOpt::OptionalValue, 
+  options.enroll("cca-path", GetLongOpt::OptionalValue,
                  "cca component path", "");
   options.enroll("cca-load", GetLongOpt::OptionalValue,
@@ -261 +229 @@
   int optind = options.parse(argc, argv);
 
-  const char *output = options.retrieve("o");
-  ostream *outstream = 0;
+  return optind;
+}
+
+/** Checks for each known option and acts accordingly.
+ *
+ * \param options option structure
+ * \param *output name of outputfile on return
+ * \param *outstream open output stream on return
+ */
+void ComputeOptions(GetLongOpt &options, const char *&output, ostream *&outstream)
+{
+  output = options.retrieve("o");
+  outstream = 0;
   if (output != 0) {
     outstream = new ofstream(output);
@@ -276 +255 @@
     exit(0);
   }
-  
+
   if (options.retrieve("v")) {
     ExEnv::out0()
@@ -284 +263 @@
     exit(0);
   }
-  
+
   if (options.retrieve("w")) {
     ExEnv::out0()
@@ -293 +272 @@
     exit(0);
   }
-  
+
   if (options.retrieve("L")) {
     ExEnv::out0()
@@ -305 +284 @@
   // set the working dir
   if (strcmp(options.retrieve("W"),"."))
-    chdir(options.retrieve("W"));
+    int retval = chdir(options.retrieve("W"));
 
   // check that n and f/o are not given at the same time
@@ -311 +290 @@
     throw invalid_argument("-n must not be given with -f or -o");
   }
+}
+
+int
+try_main(int argc, char *argv[])
+{
+  //trash_stack();
+
+  KeyValValueboolean truevalue(1), falsevalue(0);
+  int i;
+  const char *devnull = "/dev/null";
+  atexit(clean_up);
+
+#ifdef HAVE_FEENABLEEXCEPT
+  // this uses a glibc extension to trap on individual exceptions
+# ifdef FE_DIVBYZERO
+  feenableexcept(FE_DIVBYZERO);
+# endif
+# ifdef FE_INVALID
+  feenableexcept(FE_INVALID);
+# endif
+# ifdef FE_OVERFLOW
+  feenableexcept(FE_OVERFLOW);
+# endif
+#endif
+
+#ifdef HAVE_FEDISABLEEXCEPT
+  // this uses a glibc extension to not trap on individual exceptions
+# ifdef FE_UNDERFLOW
+  fedisableexcept(FE_UNDERFLOW);
+# endif
+# ifdef FE_INEXACT
+  fedisableexcept(FE_INEXACT);
+# endif
+#endif
+
+#if defined(HAVE_SETRLIMIT)
+  struct rlimit rlim;
+  rlim.rlim_cur = 0;
+  rlim.rlim_max = 0;
+  setrlimit(RLIMIT_CORE,&rlim);
+#endif
+
+  ExEnv::init(argc, argv);
+
+  // parse commandline options
+  GetLongOpt options;
+  int optind = ParseOptions(options, argc, argv);
+  const char *output = 0;
+  ostream *outstream = 0;
+  ComputeOptions(options, output, outstream);
 
   // initialize keyval input
@@ -327 +356 @@
 
   // get the message group.  first try the commandline and environment
+  Ref<MessageGrp> grp;
+#if defined(HAVE_MPI) && defined(ALWAYS_USE_MPI)
+  grp = new MPIMessageGrp(&argc, &argv);
+#endif
   if (grp.null()) grp = MessageGrp::initial_messagegrp(argc, argv);
   if (grp.nonnull())
@@ -962 +995 @@
      ExEnv::out0() << "The AO density matrix is ";
      wfn->ao_density()->print(ExEnv::out0());
+     ExEnv::out0() << "The natural density matrix is ";
+     wfn->natural_density()->print(ExEnv::out0());
      ExEnv::out0() << "The Gaussian basis is " << wfn->basis()->name() << endl;
      ExEnv::out0() << "The Gaussians sit at the following centers: " << endl;
@@ -970 +1005 @@
        ExEnv::out0() << endl;
      }
-     // GaussianShell is the actual orbital functionas it seems ...
-     ExEnv::out0() << "There are the following Gaussian Shells: " << endl;
-     for (int nr = 0; nr< wfn->basis()->nshell(); ++nr) {
-       ExEnv::out0() << "Shell nr. " << nr << ": ";
-       (*wfn->basis())[nr].print(ExEnv::out0());
+     // GaussianShell is the actual orbital functions it seems ...
+     //ExEnv::out0() << "There are the following Gaussian Shells: " << endl;
+     SCVector3 r;
+     r.x() = r.y() = r.z() = 10;
+     ExEnv::out0() << "We get the following value at " << r << "." << endl;
+     Ref<Integral> intgrl = Integral::get_default_integral();
+     GaussianBasisSet::ValueData vdat(wfn->basis(), integral);
+     ExEnv::out0() << "Value at (10,10,10) is " << EvaluateDensity(r, intgrl, vdat, wfn) << endl;
+     boost::function<double (SCVector3 &r)> evaluator =
+         boost::bind(&EvaluateDensity, _1, boost::ref(intgrl), boost::ref(vdat), boost::ref(wfn));
+     ExEnv::out0() << "Check against values at " << r << "." << endl;
+     int nbasis = wfn->basis()->nbasis();
+     double *b_val = new double[nbasis];
+     wfn->basis()->values(r, &vdat, b_val);
+     for (int i=0; i<nbasis; i++) {
+       //ExEnv::out0() << "Shell nr. " << nr << ": ";
+       ExEnv::out0() << "Value at (10,10,10) is " << b_val[i] << endl;
      }
+     // perform test integration of density
+     double delta = 1.;
+     double sum = 0.;
+     for (r.x() = -10. ; r.x() < 10.; r.x() += delta)
+       for (r.y() = -10. ; r.y() < 10.; r.y() += delta)
+         for (r.z() = -10. ; r.z() < 10.; r.z() += delta) {
+           wfn->basis()->values(r, &vdat, b_val);
+           for (int i=0; i<nbasis; i++)
+             sum += wfn->ao_density()->get_element(i,i)*b_val[i];
+         }
+     sum /= pow(20/delta,3);
+     ExEnv::out0() << "Sum over domain [0:20]^3 with " << delta << " delta is " << sum << "." << endl;
+     delete[] b_val;
   }
 
@@ -1016 +1076 @@
 }
 
+double EvaluateDensity(SCVector3 &r, Ref<Integral> &intgrl, GaussianBasisSet::ValueData &vdat, Ref<Wavefunction> &wfn)
+{
+  ExEnv::out0() << "We get the following values at " << r << "." << endl;
+  int nbasis = wfn->basis()->nbasis();
+  double *b_val = new double[nbasis];
+  wfn->basis()->values(r, &vdat, b_val);
+  double sum=0.;
+  for (int i=0; i<nbasis; i++)
+    sum += b_val[i];
+  delete[] b_val;
+  return sum;
+}
+
 int
 main(int argc, char *argv[])