Index: ThirdParty/Makefile.am
===================================================================
--- ThirdParty/Makefile.am	(revision ba1152f3b83c0c86400742bd10948d2a5729110b)
+++ ThirdParty/Makefile.am	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -14,4 +14,4 @@
 SUBDIRS += JobMarket
 endif
-SUBDIRS += LinearAlgebra
+SUBDIRS += LinearAlgebra mpqc_open
 
Index: ThirdParty/mpqc_open/src/bin/mpqc/Makefile.am
===================================================================
--- ThirdParty/mpqc_open/src/bin/mpqc/Makefile.am	(revision ba1152f3b83c0c86400742bd10948d2a5729110b)
+++ ThirdParty/mpqc_open/src/bin/mpqc/Makefile.am	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -1,34 +1,81 @@
-AM_CPPFLAGS = -I$(top_srcdir)/include -I$(top_srcdir)/src/lib
+AM_CPPFLAGS = \
+	-I$(top_srcdir)/include -I$(top_srcdir)/src/lib \
+	-DHAVE_JOBMARKET -I$(top_srcdir)/../JobMarket/src -I$(top_builddir)/../JobMarket -I$(top_srcdir)/../../src -I$(top_srcdir)/../LinearAlgebra/src -I$(top_srcdir)/../CodePatterns/src $(BOOST_SYSTEM_CFLAGS) \
+	-DHAVE_MPQCDATA -I$(top_builddir)/../../src
 
 bin_PROGRAMS = molecuilder_mpqc
 
-noinst_LTLIBRARIES = libmpqc.la
+lib_LTLIBRARIES = \
+	libmolecuilder_mpqc.la
 
-libmpqc_la_SOURCES = \
+noinst_LTLIBRARIES = \
+	libmolecuilder_mpqc_extract_dummy.la
+
+LIBMPQCSOURCES = \
+	mpqc.cc \
 	mpqcin.cc \
+	parse.cc \
+	scan.cc
+
+LIBMPQCHEADERS = \
+	mpqc.h \
 	mpqcin.h \
-	parse.cc \
 	parse.h \
-	scan.cc \
 	scdirlist.h
 
+libmolecuilder_mpqc_extract_dummy_la_SOURCES = \
+	mpqc_extract_dummy.cc \
+	mpqc_extract.h
+
+libmolecuilder_mpqc_la_includedir = $(includedir)/molecuilder_mpqc/
+libmolecuilder_mpqc_la_CPPFLAGS = $(AM_CPPFLAGS)
+libmolecuilder_mpqc_la_LDFLAGS = $(AM_LDFLAGS)
+	
+libmolecuilder_mpqc_la_LIBADD = \
+	../../lib/libSCmbpt.la ../../lib/libSCscf.la ../../lib/libSCdft.la ../../lib/libSCwfn.la ../../lib/libSCpsi.la ../../lib/libSCsolvent.la ../../lib/libSCintv3.la ../../lib/libSCbasis.la ../../lib/libSCoint3.la ../../lib/libSCmolecule.la ../../lib/libSCisosurf.la ../../lib/libSCoptimize.la ../../lib/libSCsymmetry.la ../../lib/libSCscmat.la ../../lib/libSCrender.la ../../lib/libSCgroup.la ../../lib/libSCmisc.la ../../lib/libSCstate.la ../../lib/libSCkeyval.la ../../lib/libSCclass.la ../../lib/libSCcontainer.la ../../lib/libSCref.la ../../lib/libSCoptions.la
+
+nobase_libmolecuilder_mpqc_la_include_HEADERS = ${LIBMPQCHEADERS}
+
+## Define the source file list for the "libexample-@MOLECUILDER_API_VERSION@.la"
+## target.  Note that @MOLECUILDER_API_VERSION@ is not interpreted by Automake and
+## will therefore be treated as if it were literally part of the target name,
+## and the variable name derived from that.
+## The file extension .cc is recognized by Automake, and makes it produce
+## rules which invoke the C++ compiler to produce a libtool object file (.lo)
+## from each source file.  Note that it is not necessary to list header files
+## which are already listed elsewhere in a _HEADERS variable assignment.
+libmolecuilder_mpqc_la_SOURCES = ${LIBMPQCSOURCES}
+
+## Instruct libtool to include ABI version information in the generated shared
+## library file (.so).  The library ABI version is defined in configure.ac, so
+## that all version information is kept in one place.
+libmolecuilder_mpqc_la_LDFLAGS += -version-info $(MPQC_SO_VERSION)
+
+## The generated configuration header is installed in its own subdirectory of
+## $(libdir).  The reason for this is that the configuration information put
+## into this header file describes the target platform the installed library
+## has been built for.  Thus the file must not be installed into a location
+## intended for architecture-independent files, as defined by the Filesystem
+## Hierarchy Standard (FHS).
+## The nodist_ prefix instructs Automake to not generate rules for including
+## the listed files in the distribution on 'make dist'.  Files that are listed
+## in _HEADERS variables are normally included in the distribution, but the
+## configuration header file is generated at configure time and should not be
+## shipped with the source tarball.
+libmolecuilder_mpqc_la_libincludedir = $(libdir)/molecuilder_mpqc/include
+nodist_libmolecuilder_mpqc_la_libinclude_HEADERS = $(top_builddir)/src/lib/scconfig.h
+
+## Install the generated pkg-config file (.pc) into the expected location for
+## architecture-dependent package configuration information.  Occasionally,
+## pkg-config files are also used for architecture-independent data packages,
+## in which case the correct install location would be $(datadir)/pkgconfig.
+#pkgconfigdir = $(libdir)/pkgconfig
+#pkgconfig_DATA = $(top_builddir)/molecuilder_mpqc.pc
+
 molecuilder_mpqc_SOURCES = \
-	mpqc.cc
+	mpqc_main.cc
 
-molecuilder_mpqc_CPPFLAGS = \
-	$(AM_CPPFLAGS) \
-	-DHAVE_JOBMARKET -I$(top_srcdir)/../JobMarket/src -I$(top_builddir)/../JobMarket -I$(top_srcdir)/../../src -I$(top_srcdir)/../LinearAlgebra/src $(BOOST_SYSTEM_CFLAGS) \
-	-DHAVE_MPQCDATA -I$(top_builddir)/../../src
-
-molecuilder_mpqc_LDFLAGS = \
-	$(AM_LDFLAGS) \
-	-L$(top_builddir)/../JobMarket/src/JobMarket/.libs -Wl,-rpath,$(top_builddir)/../JobMarket/src/JobMarket/.libs \
-	-L$(top_builddir)/../../src/.libs -Wl,-rpath,$(top_builddir)/../../src/.libs \
-	$(BOOST_SYSTEM_LDFLAGS)
 molecuilder_mpqc_LDADD = \
-	libmpqc.la ../../lib/libSCmbpt.la ../../lib/libSCscf.la ../../lib/libSCdft.la ../../lib/libSCwfn.la ../../lib/libSCpsi.la ../../lib/libSCsolvent.la ../../lib/libSCintv3.la ../../lib/libSCbasis.la ../../lib/libSCoint3.la ../../lib/libSCmolecule.la ../../lib/libSCisosurf.la ../../lib/libSCoptimize.la ../../lib/libSCsymmetry.la ../../lib/libSCscmat.la ../../lib/libSCrender.la ../../lib/libSCgroup.la ../../lib/libSCmisc.la ../../lib/libSCstate.la ../../lib/libSCkeyval.la ../../lib/libSCclass.la ../../lib/libSCcontainer.la ../../lib/libSCref.la ../../lib/libSCoptions.la \
-	-lJobMarketPoolWorker -lJobMarket  \
-	-lMolecuilderJobs -lMolecuilderFragmentationSummation \
-	-lboost_serialization
+	libmolecuilder_mpqc.la libmolecuilder_mpqc_extract_dummy.la
 if COND_LIBINT
 molecuilder_mpqc_LDADD += -lint
Index: ThirdParty/mpqc_open/src/bin/mpqc/mpqc.cc
===================================================================
--- ThirdParty/mpqc_open/src/bin/mpqc/mpqc.cc	(revision ba1152f3b83c0c86400742bd10948d2a5729110b)
+++ ThirdParty/mpqc_open/src/bin/mpqc/mpqc.cc	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -123,33 +123,10 @@
 #endif
 
-#ifdef HAVE_JOBMARKET
-// include headers that implement a archive in simple text format
-// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
-#include <boost/archive/text_oarchive.hpp>
-#include <boost/archive/text_iarchive.hpp>
-#include <boost/bind.hpp>
-#include <boost/tokenizer.hpp>
-
-#include "JobMarket/Results/FragmentResult.hpp"
-#include "JobMarket/poolworker_main.hpp"
-
-#include "chemistry/qc/scf/scfops.h"
-
-#ifdef HAVE_MPQCDATA
-#include "Jobs/MPQCJob.hpp"
-#include "Fragmentation/Summation/Containers/MPQCData.hpp"
-
-#include <chemistry/qc/basis/obint.h>
-#include <chemistry/qc/basis/symmint.h>
-#endif
-
-#include <algorithm>
-#include <stdlib.h>
-#endif
-
 using namespace std;
 using namespace sc;
 
 #include "mpqcin.h"
+#include "mpqc.h"
+#include "mpqc_extract.h"
 
 //////////////////////////////////////////////////////////////////////////
@@ -178,13 +155,17 @@
 }
 
-static void
-clean_up(void)
-{
-  MemoryGrp::set_default_memorygrp(0);
-  ThreadGrp::set_default_threadgrp(0);
-  SCMatrixKit::set_default_matrixkit(0);
-  Integral::set_default_integral(0);
-  RegionTimer::set_default_regiontimer(0);
-}
+namespace detail {
+
+  void
+  clean_up(void)
+  {
+    ::MemoryGrp::set_default_memorygrp(0);
+    ::ThreadGrp::set_default_threadgrp(0);
+    ::SCMatrixKit::set_default_matrixkit(0);
+    ::Integral::set_default_integral(0);
+    ::RegionTimer::set_default_regiontimer(0);
+  }
+
+} /* namespace detail */
 
 static void
@@ -335,25 +316,4 @@
 }
 
-/** Temporary structure for storing information from command-line
- *
- * This structure has been introduced to gather the various calls to GetLongOpts
- * at one (initial) place and to abstract it from the source of command-lines.
- * This temporary object can be set by other means, too. I.e. we become
- * independent of usage in command-line programs.
- */
-struct OptionValues {
-  const char *keyvalue; // option "k"
-  const char *debug; // option ""
-  int limit; // option "l"
-  const char *check; // option "c"
-  const char *simple_input; // option "i"
-  string executablename;
-
-#ifdef HAVE_CHEMISTRY_CCA
-  string cca_load;  // option "cca-path"
-  string cc_path; // option "cca-load"
-#endif
-};
-
 /** Parse remainder options not treated by ComputeOptions() into temporary storage.
  *
@@ -364,5 +324,5 @@
 void parseRemainderOptions(
     GetLongOpt &options,
-    struct OptionValues &values,
+    detail::OptionValues &values,
     int argc,
     char **argv)
@@ -655,5 +615,5 @@
 /** Parse the input configuration from char array into keyvalue container.
  *
- * \param parsedkv key value container to foll
+ * \param parsedkv key value container to fill
  * \param values temporary options value structure
  * \param in_char_array char array with input file
@@ -663,5 +623,5 @@
 void parseIntoKeyValue(
     Ref<ParsedKeyVal> &parsedkv,
-    struct OptionValues &values,
+    detail::OptionValues &values,
     char *&in_char_array,
     int use_simple_input)
@@ -699,5 +659,5 @@
     Ref<MessageGrp> &grp,
     Ref<ParsedKeyVal> &parsedkv,
-    struct OptionValues &values,
+    detail::OptionValues &values,
     const char *&input,
     const char *&generic_input,
@@ -804,5 +764,5 @@
 void prepareCCA(
     Ref<KeyVal> &keyval,
-    struct OptionValues &values
+    detail::OptionValues &values
     )
 {
@@ -851,5 +811,5 @@
     Ref<MessageGrp> &grp,
     Ref<Debugger> &debugger,
-    struct OptionValues &values)
+    detail::OptionValues &values)
 {
   debugger << keyval->describedclassvalue("debug");
@@ -971,5 +931,5 @@
 int checkBasisSetLimit(
     Ref<MolecularEnergy> &mole,
-    struct OptionValues &values
+    detail::OptionValues &values
     )
 {
@@ -1182,17 +1142,4 @@
 }
 
-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;
-}
-
 void init()
 {
@@ -1200,5 +1147,5 @@
 
   int i;
-  atexit(clean_up);
+  atexit(detail::clean_up);
 
 #ifdef HAVE_FEENABLEEXCEPT
@@ -1233,69 +1180,50 @@
 }
 
-/** Finds the region index to a given timer region name.
- *
- * @param nregion number of regions
- * @param region_names array with name of each region
- * @param name name of desired region
- * @return index of desired region in array
- */
-int findTimerRegion(const int &nregion, const char **&region_names, const char *name)
-{
-  size_t region=0;
-  for (;region<nregion;++region) {
-    //std::cout << "Comparing " << region_names[region] << " and " << name << "." << std::endl;
-    if (strcmp(region_names[region], name) == 0)
-      break;
-  }
-  if (region == nregion)
-    region = 0;
-  return region;
-}
-
 /** Performs the main work to calculate the ground state energies, gradients, etc.
  *
- * @param grp message group
- * @param values temporary value storage for parsed command-line
- * @param input input file name
- * @param generic_input generic input file name
- * @param in_char_array either NULL or contains char array with read input
+ * @param _initvalues struct with all state variables
  * @param argc argument count
  * @param argv argument array
- */
-void mainFunction(
-    Ref<MessageGrp> grp,
-    struct OptionValues &values,
-    const char *&output,
-    const char *&input,
-    const char *&generic_input,
-    char *&in_char_array,
+ * @param data void ptr to extract structure
+ */
+void mpqc::mainFunction(
+    mpqc::InitValues &_initvalues,
     int argc,
-    char **argv
-#ifdef HAVE_MPQCDATA
-    , MPQCData &data
-#endif
+    char **argv,
+    void *data
     )
 {
   // get the basename for output files
-  setOutputBaseName(input, output);
+  setOutputBaseName(_initvalues.input, _initvalues.output);
 
   // parse input into keyvalue container
   Ref<ParsedKeyVal> parsedkv;
   int use_simple_input = 0; // default is object-oriented if in_char_array is given
-  if (!in_char_array) // obtain from file
-    parseInputfile(grp, parsedkv, values, input, generic_input, in_char_array, use_simple_input);
-  parseIntoKeyValue(parsedkv, values, in_char_array, use_simple_input);
-  delete[] in_char_array;
+  if (!_initvalues.in_char_array) // obtain from file
+    parseInputfile(
+        _initvalues.grp,
+        parsedkv,
+        _initvalues.values,
+        _initvalues.input,
+        _initvalues.generic_input,
+        _initvalues.in_char_array,
+        use_simple_input);
+  parseIntoKeyValue(
+      parsedkv,
+      _initvalues.values,
+      _initvalues.in_char_array,
+      use_simple_input);
+//  delete[] in_char_array;
 
   // prefix parsed values wit "mpqc"
-  if (values.keyvalue) parsedkv->verbose(1);
+  if (_initvalues.values.keyvalue) parsedkv->verbose(1);
   Ref<KeyVal> keyval = new PrefixKeyVal(parsedkv.pointer(),"mpqc");
 
   // set up output classes
-  setupSCFormIO(grp);
+  setupSCFormIO(_initvalues.grp);
 
   // initialize timing for mpqc
   Ref<RegionTimer> tim;
-  initTimings(grp, keyval, tim);
+  initTimings(_initvalues.grp, keyval, tim);
   
   // announce ourselves
@@ -1311,6 +1239,6 @@
 
   ExEnv::out0() << indent
-       << "Using " << grp->class_name()
-       << " for message passing (number of nodes = " << grp->n() << ")." << endl
+       << "Using " << _initvalues.grp->class_name()
+       << " for message passing (number of nodes = " << _initvalues.grp->n() << ")." << endl
        << indent
        << "Using " << thread->class_name()
@@ -1320,12 +1248,12 @@
        << " for distributed shared memory." << endl
        << indent
-       << "Total number of processors = " << grp->n() * thread->nthread() << endl;
+       << "Total number of processors = " << _initvalues.grp->n() * thread->nthread() << endl;
 
   // prepare CCA if available
-  prepareCCA(keyval, values);
+  prepareCCA(keyval, _initvalues.values);
 
   // now set up the debugger
   Ref<Debugger> debugger;
-  setupDebugger(keyval, grp, debugger, values);
+  setupDebugger(keyval, _initvalues.grp, debugger, _initvalues.values);
 
   // now check to see what matrix kit to use
@@ -1370,8 +1298,8 @@
 
   // read in restart file if we do restart
-  performRestart(keyval, grp, opt, mole, restartfile);
+  performRestart(keyval, _initvalues.grp, opt, mole, restartfile);
 
   // setup molecule checkpoint file
-  setMolecularCheckpointFile(keyval, grp, mole, mole_ckpt_file);
+  setMolecularCheckpointFile(keyval, _initvalues.grp, mole, mole_ckpt_file);
   delete[] mole_ckpt_file;
 
@@ -1379,5 +1307,5 @@
   if (checkpoint && opt.nonnull()) {
     opt->set_checkpoint();
-    if (grp->me() == 0) opt->set_checkpoint_file(ckptfile);
+    if (_initvalues.grp->me() == 0) opt->set_checkpoint_file(ckptfile);
     else opt->set_checkpoint_file(devnull);
   }
@@ -1390,5 +1318,5 @@
 
   // check basis set limit
-  const int check = checkBasisSetLimit(mole, values);
+  const int check = checkBasisSetLimit(mole, _initvalues.values);
   if (check) {
     ExEnv::out0() << endl << indent
@@ -1426,5 +1354,5 @@
   // requiring optimized geometries, etc.
   if (renderer.null()) {
-    if (parsedkv.nonnull()) print_unseen(parsedkv, input);
+    if (parsedkv.nonnull()) print_unseen(parsedkv, _initvalues.input);
     keyval = 0;
     parsedkv = 0;
@@ -1465,5 +1393,5 @@
 
   // save this before doing the frequency stuff since that obsoletes the
-  saveState(wfn_file, savestate, opt, grp, mole, molname, ckptfile);
+  saveState(wfn_file, savestate, opt, _initvalues.grp, mole, molname, ckptfile);
   
   // Frequency calculation.
@@ -1473,5 +1401,5 @@
 
   if (renderer.nonnull()) {
-    renderObjects(renderer, keyval, tim, grp);
+    renderObjects(renderer, keyval, tim, _initvalues.grp);
 
     Ref<MolFreqAnimate> molfreqanim;
@@ -1489,5 +1417,5 @@
       mole->print(ExEnv::out0());
 
-    saveToPdb(do_pdb, grp, mole, molname);
+    saveToPdb(do_pdb, _initvalues.grp, mole, molname);
 
   }
@@ -1497,492 +1425,16 @@
          << endl;
   }
-  if (parsedkv.nonnull()) print_unseen(parsedkv, input);
+  if (parsedkv.nonnull()) print_unseen(parsedkv, _initvalues.input);
 
   // here, we may gather the results
   // we start to fill the MPQC_Data object
-  if (tim.nonnull()) tim->enter("gather");
-  {
-     Ref<Wavefunction> wfn;
-     wfn << mole;
-//     ExEnv::out0() << "The number of atomic orbitals: " << wfn->ao_dimension()->n() << endl;
-//     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;
-//     for (int nr = 0; nr< wfn->basis()->ncenter(); ++nr) {
-//       ExEnv::out0() << nr << " basis function has its center at ";
-//       for (int i=0; i < 3; ++i)
-//           ExEnv::out0() << wfn->basis()->r(nr,i) << "\t";
-//       ExEnv::out0() << endl;
-//     }
-     // store accuracies
-     data.accuracy = mole->value_result().actual_accuracy();
-     data.desired_accuracy = mole->value_result().desired_accuracy();
-     // print the energy
-     data.energies.total = wfn->energy();
-     data.energies.nuclear_repulsion = wfn->nuclear_repulsion_energy();
-     {
-       CLHF *clhf = dynamic_cast<CLHF*>(wfn.pointer());
-       if (clhf != NULL) {
-         double ex, ec;
-         clhf->two_body_energy(ec, ex);
-         data.energies.electron_coulomb = ec;
-         data.energies.electron_exchange = ex;
-         clhf = NULL;
-       } else {
-         ExEnv::out0() << "INFO: There is no direct CLHF information available." << endl;
-         data.energies.electron_coulomb = 0.;
-         data.energies.electron_exchange = 0.;
-       }
-     }
-     SCF *scf = NULL;
-     {
-       MBPT2 *mbpt2 = dynamic_cast<MBPT2*>(wfn.pointer());
-       if (mbpt2 != NULL) {
-         data.energies.correlation = mbpt2->corr_energy();
-         scf = mbpt2->ref().pointer();
-         CLHF *clhf = dynamic_cast<CLHF*>(scf);
-         if (clhf != NULL) {
-           double ex, ec;
-           clhf->two_body_energy(ec, ex);
-           data.energies.electron_coulomb = ec;
-           data.energies.electron_exchange = ex;
-           clhf = NULL;
-         } else {
-           ExEnv::out0() << "INFO: There is no reference CLHF information available either." << endl;
-           data.energies.electron_coulomb = 0.;
-           data.energies.electron_exchange = 0.;
-         }
-         mbpt2 = 0;
-       } else {
-         ExEnv::out0() << "INFO: There is no MBPT2 information available." << endl;
-         data.energies.correlation = 0.;
-         scf = dynamic_cast<SCF*>(wfn.pointer());
-         if (scf == NULL)
-           abort();
-       }
-     }
-     {
-       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
-
-       RefSymmSCMatrix t = scf->overlap();
-       RefSymmSCMatrix cl_dens_ = scf->ao_density();
-
-       SCFEnergy *eop = new SCFEnergy;
-       eop->reference();
-       if (t.dim()->equiv(cl_dens_.dim())) {
-         Ref<SCElementOp2> op = eop;
-         t.element_op(op,cl_dens_);
-         op=0;
-       }
-       eop->dereference();
-
-       data.energies.overlap = eop->result();
-
-       delete eop;
-       t = 0;
-       cl_dens_ = 0;
-     }
-     {
-       // taken from Wavefunction::core_hamiltonian()
-       RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
-       hao.assign(0.0);
-       Ref<PetiteList> pl = scf->integral()->petite_list();
-       Ref<SCElementOp> hc =
-         new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->kinetic(), pl));
-       hao.element_op(hc);
-       hc=0;
-
-       RefSymmSCMatrix h(scf->so_dimension(), scf->basis_matrixkit());
-       pl->symmetrize(hao,h);
-
-       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
-       RefSymmSCMatrix cl_dens_ = scf->ao_density();
-
-       SCFEnergy *eop = new SCFEnergy;
-       eop->reference();
-       if (h.dim()->equiv(cl_dens_.dim())) {
-         Ref<SCElementOp2> op = eop;
-         h.element_op(op,cl_dens_);
-         op=0;
-       }
-       eop->dereference();
-
-       data.energies.kinetic = 2.*eop->result();
-
-       delete eop;
-       hao = 0;
-       h = 0;
-       cl_dens_ = 0;
-     }
-     {
-       // set to potential energy between nuclei and electron charge distribution
-       RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
-       hao.assign(0.0);
-       Ref<PetiteList> pl = scf->integral()->petite_list();
-       Ref<SCElementOp> hc =
-         new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->nuclear(), pl));
-       hao.element_op(hc);
-       hc=0;
-
-       RefSymmSCMatrix h(scf->so_dimension(), scf->basis_matrixkit());
-       pl->symmetrize(hao,h);
-
-       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
-       RefSymmSCMatrix cl_dens_ = scf->ao_density();
-
-       SCFEnergy *eop = new SCFEnergy;
-       eop->reference();
-       if (h.dim()->equiv(cl_dens_.dim())) {
-         Ref<SCElementOp2> op = eop;
-         h.element_op(op,cl_dens_);
-         op=0;
-       }
-       eop->dereference();
-
-       data.energies.hcore = 2.*eop->result();
-
-       delete eop;
-       hao = 0;
-       h = 0;
-       cl_dens_ = 0;
-     }
-     ExEnv::out0() << "total is " << data.energies.total << endl;
-     ExEnv::out0() << "nuclear_repulsion is " << data.energies.nuclear_repulsion << endl;
-     ExEnv::out0() << "electron_coulomb is " << data.energies.electron_coulomb << endl;
-     ExEnv::out0() << "electron_exchange is " << data.energies.electron_exchange << endl;
-     ExEnv::out0() << "correlation is " << data.energies.correlation << endl;
-     ExEnv::out0() << "overlap is " << data.energies.overlap << endl;
-     ExEnv::out0() << "kinetic is " << data.energies.kinetic << endl;
-     ExEnv::out0() << "hcore is " << data.energies.hcore << endl;
-     ExEnv::out0() << "sum is " <<
-         data.energies.nuclear_repulsion
-         + data.energies.electron_coulomb
-         + data.energies.electron_exchange
-         + data.energies.correlation
-         + data.energies.kinetic
-         + data.energies.hcore
-         << endl;
-
-     ExEnv::out0() << endl << indent
-          << scprintf("Value of the MolecularEnergy: %15.10f",
-                      mole->energy())
-          << endl;
-     // print the gradient
-     RefSCVector grad;
-     if (mole->gradient_result().computed()) {
-       grad = mole->gradient_result().result_noupdate();
-     }
-     // gradient calculation needs to be activated in the configuration
-     // some methods such as open shell MBPT2 do not allow for gradient calc.
-//     else {
-//       grad = mole->gradient();
-//     }
-     if (grad.nonnull()) {
-       data.forces.resize(grad.dim()/3);
-       for (int j=0;j<grad.dim()/3; ++j) {
-         data.forces[j].resize(3, 0.);
-       }
-       ExEnv::out0() << "Gradient of the MolecularEnergy:" << std::endl;
-       for (int j=0;j<grad.dim()/3; ++j) {
-         ExEnv::out0() << "\t";
-         for (int i=0; i< 3; ++i) {
-           data.forces[j][i] = grad[3*j+i];
-           ExEnv::out0() << grad[3*j+i] << "\t";
-         }
-         ExEnv::out0() << endl;
-       }
-     } else {
-       ExEnv::out0() << "INFO: There is no gradient information available." << endl;
-     }
-     grad = NULL;
-
-     {
-       // eigenvalues (this only works if we have a OneBodyWavefunction, i.e. SCF procedure)
-       // eigenvalues seem to be invalid for unrestricted SCF calculation 
-       // (see UnrestrictedSCF::eigenvalues() implementation)
-       UnrestrictedSCF *uscf = dynamic_cast<UnrestrictedSCF*>(wfn.pointer());
-       if ((scf != NULL) && (uscf == NULL)) {
-         const double scfernergy = scf->energy();
-         RefDiagSCMatrix evals = scf->eigenvalues();
-
-         ExEnv::out0() << "Eigenvalues:" << endl;
-         for(int i=0;i<wfn->oso_dimension(); ++i) {
-           data.energies.eigenvalues.push_back(evals(i));
-           ExEnv::out0() << i << "th eigenvalue is " << evals(i) << endl;
-         }
-       } else {
-           ExEnv::out0() << "INFO: There is no eigenvalue information available." << endl;
-       }
-     }
-     // we do sample the density only on request
-       {
-         // fill positions and charges (NO LONGER converting from bohr radii to angstroem)
-         const double AtomicLengthToAngstroem = 1.;//0.52917721;
-         data.positions.reserve(wfn->molecule()->natom());
-         data.atomicnumbers.reserve(wfn->molecule()->natom());
-         data.charges.reserve(wfn->molecule()->natom());
-         for (int iatom=0;iatom < wfn->molecule()->natom(); ++iatom) {
-           data.atomicnumbers.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)
-             pos[j] = wfn->molecule()->r(iatom, j)*AtomicLengthToAngstroem;
-           data.positions.push_back(pos);
-         }
-         ExEnv::out0() << "We have "
-             << data.positions.size() << " positions and "
-             << data.charges.size() << " charges." << endl;
-       }
-       if (data.DoLongrange) {
-       if (data.sampled_grid.level != 0)
-       {
-         // we now need to sample the density on the grid
-         // 1. get max and min over all basis function positions
-         assert( scf->basis()->ncenter() > 0 );
-         SCVector3 bmin( scf->basis()->r(0,0), scf->basis()->r(0,1), scf->basis()->r(0,2) );
-         SCVector3 bmax( scf->basis()->r(0,0), scf->basis()->r(0,1), scf->basis()->r(0,2) );
-         for (int nr = 1; nr< scf->basis()->ncenter(); ++nr) {
-           for (int i=0; i < 3; ++i) {
-             if (scf->basis()->r(nr,i) < bmin(i))
-               bmin(i) = scf->basis()->r(nr,i);
-             if (scf->basis()->r(nr,i) > bmax(i))
-               bmax(i) = scf->basis()->r(nr,i);
-           }
-         }
-         ExEnv::out0() << "Basis min is at " << bmin << " and max is at " << bmax << endl;
-
-         // 2. choose an appropriately large grid
-         // we have to pay attention to capture the right amount of the exponential decay
-         // and also to have a power of two size of the grid at best
-         SCVector3 boundaryV(5.);  // boundary extent around compact domain containing basis functions
-         bmin -= boundaryV;
-         bmax += boundaryV;
-         for (size_t i=0;i<3;++i) {
-           if (bmin(i) < data.sampled_grid.begin[i])
-             bmin(i) = data.sampled_grid.begin[i];
-           if (bmax(i) > data.sampled_grid.end[i])
-             bmax(i) = data.sampled_grid.end[i];
-         }
-         // set the non-zero window of the sampled_grid
-         data.sampled_grid.setWindow(bmin.data(), bmax.data());
-
-         // for the moment we always generate a grid of full size
-         // (NO LONGER converting grid dimensions from angstroem to bohr radii)
-         const double AtomicLengthToAngstroem = 1.;//0.52917721;
-         SCVector3 min;
-         SCVector3 max;
-         SCVector3 delta;
-         size_t samplepoints[3];
-         // due to periodic boundary conditions, we don't need gridpoints-1 here
-         // TODO: in case of open boundary conditions, we need data on the right
-         // hand side boundary as well
-         const int gridpoints = data.sampled_grid.getGridPointsPerAxis();
-         for (size_t i=0;i<3;++i) {
-           min(i) = data.sampled_grid.begin_window[i]/AtomicLengthToAngstroem;
-           max(i) = data.sampled_grid.end_window[i]/AtomicLengthToAngstroem;
-           delta(i) = data.sampled_grid.getDeltaPerAxis(i)/AtomicLengthToAngstroem;
-           samplepoints[i] = data.sampled_grid.getWindowGridPointsPerAxis(i);
-         }
-         ExEnv::out0() << "Grid starts at " << min
-             << " and ends at " << max
-             << " with a delta of " << delta
-             << " to get "
-             << samplepoints[0] << ","
-             << samplepoints[1] << ","
-             << samplepoints[2] << " samplepoints."
-             << endl;
-         assert( data.sampled_grid.sampled_grid.size() == samplepoints[0]*samplepoints[1]*samplepoints[2]);
-
-         // 3. sample the atomic density
-         const double element_volume_conversion =
-             1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
-         SCVector3 r = min;
-
-         std::set<int> valence_indices;
-         RefDiagSCMatrix evals = scf->eigenvalues();
-         if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly) {
-           // find valence orbitals
-//           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;
-           int n_electrons = scf->nelectron();
-           int n_core_electrons =  wfn->molecule()->n_core_electrons();
-           std::set<double> evals_sorted;
-           {
-             int i=0;
-             double first_positive_ev = std::numeric_limits<double>::max();
-             for(i=0;i<wfn->oso_dimension(); ++i) {
-               if (evals(i) < 0.)
-                 evals_sorted.insert(evals(i));
-               else
-                 first_positive_ev = std::min(first_positive_ev, (double)evals(i));
-             }
-             // add the first positive for the distance
-             evals_sorted.insert(first_positive_ev);
-           }
-           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 might be " << last_core_ev << std::endl;
-           ExEnv::out0()  << "First valence index is " << n_core_electrons/2 << std::endl;
-           for(int i=n_core_electrons/2;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);
-           ExEnv::out0() << "# " << i << " has ev of " << evals(i) << ", occupied by " << scf->occupation(i) << std::endl;
-         }
-         RefSymmSCMatrix d2(scf->so_dimension(), scf->basis_matrixkit());
-         d2.assign(0.0);
-         d2.accumulate_transform(vector, occ);
-
-         // taken from scf::density()
-         RefSCMatrix nos
-             = scf->integral()->petite_list()->evecs_to_AO_basis(scf->natural_orbitals());
-         RefDiagSCMatrix nd = scf->natural_density();
-         GaussianBasisSet::ValueData *valdat
-           = new GaussianBasisSet::ValueData(scf->basis(), scf->integral());
-         std::vector<double>::iterator griditer = data.sampled_grid.sampled_grid.begin();
-         const int nbasis = scf->basis()->nbasis();
-         double *bs_values = new double[nbasis];
-
-         // TODO: need to take care when we have periodic boundary conditions.
-         for (size_t x = 0; x < samplepoints[0]; ++x, r.x() += delta(0)) {
-           std::cout << "Sampling now for x=" << r.x() << std::endl;
-           for (size_t y = 0; y < samplepoints[1]; ++y, r.y() += delta(1)) {
-             for (size_t z = 0; z < samplepoints[2]; ++z, r.z() += delta(2)) {
-               scf->basis()->values(r,valdat,bs_values);
-
-               // loop over natural orbitals adding contributions to elec_density
-               double elec_density=0.0;
-               for (int i=0; i<nbasis; ++i) {
-                   double tmp = 0.0;
-                   for (int j=0; j<nbasis; ++j) {
-                       tmp += d2(j,i)*bs_values[j]*bs_values[i];
-                     }
-                   elec_density += tmp;
-                 }
-               const double dens_at_r = elec_density * element_volume_conversion;
-  //             const double dens_at_r = scf->density(r) * element_volume_conversion;
-
-  //             if (fabs(dens_at_r) > 1e-4)
-  //               std::cout << "Electron density at " << r << " is " << dens_at_r << std::endl;
-               if (griditer != data.sampled_grid.sampled_grid.end())
-                 *griditer++ = dens_at_r;
-               else
-                 std::cerr << "PAST RANGE!" << std::endl;
-             }
-             r.z() = min.z();
-           }
-           r.y() = min.y();
-         }
-         delete[] bs_values;
-         delete valdat;
-         assert( griditer == data.sampled_grid.sampled_grid.end());
-         // normalization of electron charge to equal electron number
-         {
-           double integral_value = 0.;
-           const double volume_element = pow(AtomicLengthToAngstroem,3)*delta(0)*delta(1)*delta(2);
-           for (std::vector<double>::const_iterator diter = data.sampled_grid.sampled_grid.begin();
-               diter != data.sampled_grid.sampled_grid.end(); ++diter)
-               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) || (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 " << ((data.DoValenceOnly == MPQCData::DoSampleValenceOnly) ? "n_valence_electrons" : "n_electrons")
-               << " of " << n_electrons << "." << std::endl;
-           // with normalization we also get the charge right : -1.
-           for (std::vector<double>::iterator diter = data.sampled_grid.sampled_grid.begin();
-               diter != data.sampled_grid.sampled_grid.end(); ++diter)
-               *diter *= -1.*normalization;
-         }
-       }
-     }
-     scf = 0;
-  }
-  if (tim.nonnull()) tim->exit("gather");
+  if (tim.nonnull()) tim->enter("extract");
+  extractResults(mole, data);
+  if (tim.nonnull()) tim->exit("extract");
 
   if (print_timings)
     if (tim.nonnull()) tim->print(ExEnv::out0());
 
-
-  {
-    // times obtain from key "mpqc" which should be the first
-    const int nregion = tim->nregion();
-    //std::cout << "There are " << nregion << " timed regions." << std::endl;
-    const char **region_names = new const char*[nregion];
-    tim->get_region_names(region_names);
-    // find "gather"
-    size_t gather_region = findTimerRegion(nregion, region_names, "gather");
-    size_t mpqc_region = findTimerRegion(nregion, region_names, "mpqc");
-    delete[] region_names;
-
-    // get timings
-    double *cpu_time = new double[nregion];
-    double *wall_time = new double[nregion];
-    double *flops = new double[nregion];
-    tim->get_cpu_times(cpu_time);
-    tim->get_wall_times(wall_time);
-    tim->get_flops(flops);
-    if (cpu_time != NULL) {
-      data.times.total_cputime = cpu_time[mpqc_region];
-      data.times.gather_cputime = cpu_time[gather_region];
-    }
-    if (wall_time != NULL) {
-      data.times.total_walltime = wall_time[mpqc_region];
-      data.times.gather_walltime = wall_time[gather_region];
-    }
-    if (flops != NULL) {
-      data.times.total_flops = flops[mpqc_region];
-      data.times.gather_flops = flops[gather_region];
-    }
-    delete[] cpu_time;
-    delete[] wall_time;
-    delete[] flops;
-  }
+  extractTimings(tim, data);
 
   delete[] molname;
@@ -2001,5 +1453,5 @@
   parsedkv = 0;
   memory = 0;
-  clean_up();
+  detail::clean_up();
 
 #if defined(HAVE_TIME) && defined(HAVE_CTIME)
@@ -2016,186 +1468,91 @@
 
 // static values object
-OptionValues values;
-
-#ifdef HAVE_JOBMARKET
-FragmentResult::ptr MPQCJob::Work()
-{
-  char mpqc[] = "mpqc" ;
-  char **argv = new char*[1];
-  argv[0] = &mpqc[0];
-  int argc = 1;
-//  init();
-//
-//  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);
-//  OptionValues values;
-//  parseRemainderOptions(options, values, argc, argv);
-//
-//  // get the basename for output files
-//  char filename_template[] = "mpqc_temp_XXXXXX";
-//  output = mktemp(filename_template);
-//  setOutputBaseName(NULL, output);
-
-  // now comes the actual work
-  int nfilebase = (int) inputfile.length();
-  char *in_char_array = new char[nfilebase + 1];
-  strncpy(in_char_array, inputfile.c_str(), nfilebase);
-  in_char_array[nfilebase] = '\0';
-  const char *input = 0;
-  const char *generic_input = 0;
-  Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
-  // create unique, temporary name and check whether it exists
-  const char *output = NULL;
-  char *tempfilename = NULL;
+detail::OptionValues values;
+
+
+namespace mpqc {
+
+  InitValues::InitValues() :
+      input(NULL),
+      object_input(NULL),
+      generic_input(NULL),
+      output(NULL),
+      outstream(NULL),
+      in_char_array(NULL),
+      values(::values)
+  {};
+
+  void init(InitValues &_initvalues, int argc, char *argv[]) {
+
+#if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
+    MPI_Init(&argc, &argv);
+#endif
+
+    ::init();
+
+    ExEnv::init(argc, argv);
+
+  }
+
+  void parseOptions(InitValues &_initvalues, int argc, char *argv[])
   {
-    std::ifstream test;
-    do {
-      if (output != NULL) // free buffer from last round
-        delete output;
-      char filename_template[] = "mpqc_temp_XXXXXX\0";
-      char filename_suffix[] = ".in\0";
-      tempfilename = (char *) malloc ( (strlen(filename_template)+strlen(filename_suffix)+2)*(sizeof(char)));
-      strncpy(tempfilename, mktemp(filename_template), strlen(filename_template));
-      tempfilename[strlen(filename_template)] = '\0';
-      strncat(tempfilename, filename_suffix, strlen(filename_suffix));
-      output = tempfilename;
-      //free(tempfilename); // don't free! output takes over pointer!
-      test.open(output);
-    } while (test.good()); // test whether file does not(!) exist
-    test.close();
-  }
-  // put info how to sample the density into MPQCData
-  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
-  try {
-    mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);
-  }
-  catch (SCException &e) {
-    cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
-         << e.what()
-         << endl;
-    clean_up();
-  }
-
-  //delete[] in_char_array; // is deleted in mainFunction()
-  if (output != 0) {
-    free(tempfilename);
-  }
-  delete[] argv;
-  grp = NULL;
-
-  // place into returnstream
-  std::stringstream returnstream;
-  boost::archive::text_oarchive oa(returnstream);
-  oa << data;
-
-  FragmentResult::ptr s( new FragmentResult(getId(), returnstream.str()) );
-  if (s->exitflag != 0)
-    cerr << "Job #" << s->getId() << " failed to reach desired accuracy." << endl;
-
-  return s;
-}
-#endif
-
-// we need to explicitly instantiate the serialization functions as
-// its is only serialized through its base class FragmentJob
-BOOST_CLASS_EXPORT_IMPLEMENT(MPQCJob)
+    // parse commandline options
+    int optind = ParseOptions(_initvalues.options, argc, argv);
+    ComputeOptions(_initvalues.options, _initvalues.output, _initvalues.outstream);
+    parseRemainderOptions(_initvalues.options, _initvalues.values, argc, argv);
+
+    // get input file names, either object-oriented or generic
+    getInputFileNames(_initvalues.object_input, _initvalues.generic_input, _initvalues.options, optind, argc, argv);
+    if (_initvalues.object_input)
+      _initvalues.input = _initvalues.object_input;
+    if (_initvalues.generic_input)
+      _initvalues.input = _initvalues.generic_input;
+
+    // get the message group.  first try the commandline and environment
+    getMessageGroup(_initvalues.grp, argc, argv);
+  }
+
+  void cleanup(InitValues &_initvalues) {
+    if (_initvalues.output != 0) {
+      ExEnv::set_out(&cout);
+      delete _initvalues.outstream;
+    }
+
+    _initvalues.grp = 0;
+    final_clean_up();
+
+#if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
+    // there is an MPI_Finalize() hidden somewhere else
+//    MPI_Finalize();
+#endif
+  }
+
+} /* namespace mpqc */
+
+namespace detail {
 
 int
 try_main(int argc, char *argv[])
 {
-  init();
-
-  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);
-  parseRemainderOptions(options, values, argc, argv);
-
-  // get input file names, either object-oriented or generic
-  const char *object_input = 0;
-  const char *generic_input = 0;
-  getInputFileNames(object_input, generic_input, options, optind, argc, argv);
-  const char *input = 0;
-  if (object_input) input = object_input;
-  if (generic_input) input = generic_input;
-
-  // get the message group.  first try the commandline and environment
-  Ref<MessageGrp> grp;
-  getMessageGroup(grp, argc, argv);
+  mpqc::InitValues initvalues;
+
+  mpqc::init(initvalues, argc, argv);
+
+  mpqc::parseOptions(initvalues, argc, argv);
 
   // check if we got option "-n"
   int exitflag = 0;
-#ifdef HAVE_JOBMARKET
-  if (options.retrieve("n")) {
-    /// create new argc, argv and call by splitting with tokenizer
-    std::string networkstring(options.retrieve("n"));
-    typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
-    boost::char_separator<char> sep(" ");
-    tokenizer tok(networkstring, sep);
-    // simple count because tokenizer has now size()
-    int argc_network = 0;
-    for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg)
-      ++argc_network;
-    // argv[0] is program name
-    char **argv_network = new char*[argc_network+1];
-    argv_network[0] = new char[5];
-    strcpy(argv_network[0], "mpqc");
-    // then we place each token as a new argument
-    argc_network = 1;
-    for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg){
-      const size_t strlength = (*beg).length();
-      const char *strarray = beg->c_str();
-      //std::cout << "Token " << argc_network << " is " << strarray << ", length " << strlength << endl;
-      argv_network[argc_network] = new char[strlength+1];
-      strcpy(argv_network[argc_network], strarray);
-      argv_network[argc_network][strlength] = '\0';
-      for (size_t index = 0; index < strlength; ++index)
-        if (argv_network[argc_network][index] == '+')
-          argv_network[argc_network][index] = '-';
-      ++argc_network;
-    }
-    /// and start listening for MPQCJobs
-    exitflag = poolworker_main(argc_network, argv_network);
-    /// remove the artifical [argv,argv] again
-    for (int i=0;i<argc_network;++i)
-      delete[] argv_network[i];
-    delete[] argv_network;
-  } else
-#endif
-  {
-    // if not, work on the command line input
-    char *in_char_array = 0;
-#ifdef HAVE_MPQCDATA
-    MPQCData data;
-    mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);
-#else
-    mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv);
-#endif
-  }
-
-  if (output != 0) {
-    ExEnv::set_out(&cout);
-    delete outstream;
-  }
-
-  grp = 0;
-  final_clean_up();
+  void *data = NULL;
+    mpqc::mainFunction(
+        initvalues,
+        argc, argv,
+        data);
+
+  mpqc::cleanup(initvalues);
 
   return exitflag;
 }
 
+} /* namespace detail */
 
 double EvaluateDensity(SCVector3 &r, Ref<Integral> &intgrl, GaussianBasisSet::ValueData &vdat, Ref<Wavefunction> &wfn)
@@ -2212,47 +1569,4 @@
 }
 
-int
-main(int argc, char *argv[])
-{
-#if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
-  MPI_Init(&argc, &argv);
-#endif
-  size_t exitflag = 0;
-  try {
-    exitflag = try_main(argc, argv);
-  }
-  catch (SCException &e) {
-    cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
-         << e.what()
-         << endl;
-    clean_up();
-    throw;
-  }
-  catch (bad_alloc &e) {
-    cout << argv[0] << ": ERROR: MEMORY ALLOCATION FAILED:" << endl
-         << e.what()
-         << endl;
-    clean_up();
-    throw;
-  }
-  catch (exception &e) {
-    cout << argv[0] << ": ERROR: EXCEPTION RAISED:" << endl
-         << e.what()
-         << endl;
-    clean_up();
-    throw;
-  }
-  catch (...) {
-    cout << argv[0] << ": ERROR: UNKNOWN EXCEPTION RAISED" << endl;
-    clean_up();
-    throw;
-  }
-#if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
-  // there is an MPI_Finalize() hidden somewhere else
-//  MPI_Finalize();
-#endif
-  return exitflag;
-}
-
 /////////////////////////////////////////////////////////////////////////////
 
Index: ThirdParty/mpqc_open/src/bin/mpqc/mpqc.h
===================================================================
--- ThirdParty/mpqc_open/src/bin/mpqc/mpqc.h	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
+++ ThirdParty/mpqc_open/src/bin/mpqc/mpqc.h	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -0,0 +1,101 @@
+
+#ifndef _mpqc_h
+#define _mpqc_h
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+
+#include <iosfwd>
+#include <string>
+
+#include <util/group/message.h>
+#include <util/options/GetLongOpt.h>
+
+namespace detail {
+
+  /** Temporary structure for storing information from command-line
+   *
+   * This structure has been introduced to gather the various calls to GetLongOpts
+   * at one (initial) place and to abstract it from the source of command-lines.
+   * This temporary object can be set by other means, too. I.e. we become
+   * independent of usage in command-line programs.
+   */
+  struct OptionValues {
+    const char *keyvalue; // option "k"
+    const char *debug; // option ""
+    int limit; // option "l"
+    const char *check; // option "c"
+    const char *simple_input; // option "i"
+    std::string executablename;
+
+  #ifdef HAVE_CHEMISTRY_CCA
+    string cca_load;  // option "cca-path"
+    string cc_path; // option "cca-load"
+  #endif
+  };
+
+  int try_main(int argc, char *argv[]);
+
+  void clean_up(void);
+
+} /* namespace detail */
+
+namespace mpqc {
+
+  /** This is just a helper structure to combine all variables that get
+   * initialized in init()
+   */
+  struct InitValues {
+    InitValues();
+
+    const char *input;
+    const char *object_input;
+    const char *generic_input;
+    const char *output;
+    std::ostream *outstream;
+    char *in_char_array;
+    detail::OptionValues &values;
+    sc::GetLongOpt options;
+    sc::Ref< sc::MessageGrp> grp;
+  };
+
+  /** Initializes all variables inside InitValues struct and prepares code
+   * for calculation.
+   *
+   * \param _initvalues struct holding state of variables to initialize
+   * \param argc argument count
+   * \param argv argument variables array
+   */
+  void init(InitValues &_initvalues, int argc, char *argv[]);
+
+  /** Parses command-line options, sets files and initializes message group
+   *
+   * \param _initvalues struct holding state of variables to initialize
+   * \param argc argument count
+   * \param argv argument variables array
+   */
+  void parseOptions(InitValues &_initvalues, int argc, char *argv[]);
+
+  /** Performs the main work to calculate the ground state energies, gradients, etc.
+   *
+   * @param _initvalues struct with all state variables
+   * @param argc argument count
+   * @param argv argument array
+   * @param data void ptr to extract structure
+   */
+  void mainFunction(
+      InitValues &_initvalues,
+      int argc,
+      char **argv,
+      void *data
+      );
+
+  /** Cleans up variables in initialized variables struct.
+   *
+   */
+  void cleanup(InitValues &_initvalues);
+
+} /* namespace mpqc */
+
+#endif
Index: ThirdParty/mpqc_open/src/bin/mpqc/mpqc_extract.h
===================================================================
--- ThirdParty/mpqc_open/src/bin/mpqc/mpqc_extract.h	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
+++ ThirdParty/mpqc_open/src/bin/mpqc/mpqc_extract.h	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -0,0 +1,19 @@
+#ifndef _mpqc_extract_h
+#define _mpqc_extract_h
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+
+#include <chemistry/molecule/energy.h>
+#include <util/group/pregtime.h>
+
+void extractResults(
+    sc::Ref<sc::MolecularEnergy> &mole,
+    void *_data);
+
+void extractTimings(
+    sc::Ref<sc::RegionTimer> &tim,
+    void *_data);
+
+#endif
Index: ThirdParty/mpqc_open/src/bin/mpqc/mpqc_extract_dummy.cc
===================================================================
--- ThirdParty/mpqc_open/src/bin/mpqc/mpqc_extract_dummy.cc	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
+++ ThirdParty/mpqc_open/src/bin/mpqc/mpqc_extract_dummy.cc	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -0,0 +1,44 @@
+/* Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2017 Frederik Heber. All rights reserved.
+ *
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/** \file mpqc_extract_dummy.cc
+ *
+ * Function implementations for the class vector.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+
+#include "mpqc_extract.h"
+
+using namespace sc;
+
+// dummy implementations of extractor functions
+
+void extractResults(
+    Ref<MolecularEnergy> &mole,
+    void *_data) {}
+
+void extractTimings(
+    Ref<RegionTimer> &tim,
+    void *_data) {}
Index: ThirdParty/mpqc_open/src/bin/mpqc/mpqc_main.cc
===================================================================
--- ThirdParty/mpqc_open/src/bin/mpqc/mpqc_main.cc	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
+++ ThirdParty/mpqc_open/src/bin/mpqc/mpqc_main.cc	(revision 7672284238c24f3166b6ef5766db518d837c36ad)
@@ -0,0 +1,85 @@
+//
+// mpqc_main.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of MPQC.
+//
+// MPQC is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// MPQC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with the MPQC; see the file COPYING.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+// \note This was extracted from \file mpqc.cc for refactoring into a library.
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+
+#include <iostream>
+
+#include <util/class/scexception.h>
+
+#include "mpqc.h"
+
+using namespace std;
+
+using namespace sc;
+
+int
+main(int argc, char *argv[])
+{
+  size_t exitflag = 0;
+  try {
+    exitflag = detail::try_main(argc, argv);
+  }
+  catch (SCException &e) {
+    cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
+         << e.what()
+         << endl;
+    detail::clean_up();
+    throw;
+  }
+  catch (bad_alloc &e) {
+    cout << argv[0] << ": ERROR: MEMORY ALLOCATION FAILED:" << endl
+         << e.what()
+         << endl;
+    detail::clean_up();
+    throw;
+  }
+  catch (exception &e) {
+    cout << argv[0] << ": ERROR: EXCEPTION RAISED:" << endl
+         << e.what()
+         << endl;
+    detail::clean_up();
+    throw;
+  }
+  catch (...) {
+    cout << argv[0] << ": ERROR: UNKNOWN EXCEPTION RAISED" << endl;
+    detail::clean_up();
+    throw;
+  }
+
+  return exitflag;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
