Index: src/Jobs/MPQCJob_Work.cpp
===================================================================
--- src/Jobs/MPQCJob_Work.cpp	(revision 69682af4a24e4cabdca10e5bae32c1b59a630b37)
+++ src/Jobs/MPQCJob_Work.cpp	(revision fbf005fe6f3de1a0fe7c1467b64db03e679102ae)
@@ -1,7 +1,6 @@
-/*
- * Project: MoleCuilder
+/* Project: MoleCuilder
  * Description: creates and alters molecular systems
  * Copyright (C)  2017 Frederik Heber. All rights reserved.
- * 
+ *
  *
  *   This file is part of MoleCuilder.
@@ -21,11 +20,9 @@
  */
 
-/*
- * MPQCJobi_Work.cpp
+/** \file MPCJob_Work.cpp
  *
- *  Created on: Feb 20, 2017
- *      Author: heber
+ * Function implementations for the class vector.
+ *
  */
-
 
 // include config.h
@@ -34,17 +31,99 @@
 #endif
 
+#include <sstream>
+
+#include "MPQCJob.hpp"
+
+#include "CodePatterns/Log.hpp"
+
+#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/tokenizer.hpp>
 
-#include "CodePatterns/MemDebug.hpp"
+#include "JobMarket/Results/FragmentResult.hpp"
 
-#include "MPQCJob.hpp"
+#include <fstream>
+#include <iostream>
+#include <stdlib.h>
 
-#include <sstream>
+#include <util/class/scexception.h>
+#include <util/group/message.h>
 
-#include "CodePatterns/Log.hpp"
+#include "mpqc.h"
 
+using namespace std;
+
+using namespace sc;
+
+FragmentResult::ptr MPQCJob::Work()
+{
+  char mpqc[] = "mpqc" ;
+  char **argv = new char*[1];
+  argv[0] = &mpqc[0];
+  int argc = 1;
+
+  // now comes the actual work
+  mpqc::InitValues initvalues;
+  int nfilebase = (int) inputfile.length();
+  initvalues.in_char_array = new char[nfilebase + 1];
+  strncpy(initvalues.in_char_array, inputfile.c_str(), nfilebase);
+  initvalues.in_char_array[nfilebase] = '\0';
+  initvalues.grp = MessageGrp::get_default_messagegrp();
+  // create unique, temporary name and check whether it exists
+  char *tempfilename = NULL;
+  {
+    std::ifstream test;
+    do {
+      if (initvalues.output != NULL) // free buffer from last round
+        delete initvalues.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));
+      initvalues.output = tempfilename;
+      //free(tempfilename); // don't free! output takes over pointer!
+      test.open(initvalues.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 {
+    mpqc::mainFunction(initvalues, argc, argv, static_cast<void *>(&data));
+  }
+  catch (SCException &e) {
+    cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
+         << e.what()
+         << endl;
+    detail::clean_up();
+  }
+
+  //delete[] in_char_array; // is deleted in mainFunction()
+  if (initvalues.output != 0) {
+    free(tempfilename);
+  }
+  delete[] argv;
+  initvalues.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;
+}
+#else
 FragmentResult::ptr MPQCJob::Work()
 {
@@ -67,3 +146,7 @@
   return ptr;
 }
+#endif
 
+// we need to explicitly instantiate the serialization functions as
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_IMPLEMENT(MPQCJob)
Index: src/Jobs/Makefile.am
===================================================================
--- src/Jobs/Makefile.am	(revision 69682af4a24e4cabdca10e5bae32c1b59a630b37)
+++ src/Jobs/Makefile.am	(revision fbf005fe6f3de1a0fe7c1467b64db03e679102ae)
@@ -56,23 +56,55 @@
 
 lib_LTLIBRARIES += \
-	libMolecuilderJobs.la
+	libMolecuilderJobs.la \
+	libmolecuilder_mpqc_extract.la
 noinst_LTLIBRARIES += \
 	libMolecuilderJobs_Work.la
+
+### libmolecuilder_mpqc_extract.la
+libmolecuilder_mpqc_extract_la_SOURCES = \
+	Jobs/mpqc_extract.cc \
+	$(top_srcdir)/ThirdParty/mpqc_open/src/bin/mpqc/mpqc.h
+libmolecuilder_mpqc_extract_la_LDFLAGS = \
+	$(AM_CPPFLAGS) \
+	-L$(abs_top_builddir)/../JobMarket/src/JobMarket/.libs -Wl,-rpath,$(abs_top_builddir)/../JobMarket/src/JobMarket/.libs \
+	-L$(abs_top_builddir)/../../src/.libs -Wl,-rpath,$(abs_top_builddir)/../../src/.libs \
+	$(BOOST_SYSTEM_LDFLAGS)
+libmolecuilder_mpqc_extract_la_CPPFLAGS = \
+	$(AM_CPPFLAGS) \
+	${CodePatterns_CFLAGS} ${JobMarket_CFLAGS} \
+	$(BOOST_SYSTEM_CFLAGS) \
+	-DHAVE_MPQCDATA \
+	-DHAVE_JOBMARKET \
+	-I$(top_builddir)/ThirdParty/mpqc_open/src/lib \
+	-I$(top_srcdir)/ThirdParty/mpqc_open/src/lib \
+	-I$(top_srcdir)/ThirdParty/mpqc_open/src/bin/mpqc
+libmolecuilder_mpqc_extract_la_LIBADD = \
+	libMolecuilderJobs.la \
+	libMolecuilderFragmentationSummation.la \
+	-lboost_serialization
+
+### libMolecuilderJobs_Work.la
+libMolecuilderJobs_Work_la_CPPFLAGS = \
+	$(AM_CPPFLAGS) ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} ${JobMarket_CFLAGS} -Dvmg_float=double -Dvmg_int=int $(VMG_CFLAGS) \
+	-I$(top_builddir)/ThirdParty/mpqc_open/src/lib \
+	-I$(top_srcdir)/ThirdParty/mpqc_open/src/lib \
+	-I$(top_srcdir)/ThirdParty/mpqc_open/src/bin/mpqc
+libMolecuilderJobs_Work_la_LDFLAGS = $(AM_LDFLAGS)
+if CONDJOBMARKET
+libMolecuilderJobs_Work_la_LDFLAGS += \
+  $(JobMarket_LDFLAGS)
+endif
+
+### libMolecuilderJobs.la
 libMolecuilderJobs_la_includedir = $(includedir)/MoleCuilder/
 nobase_libMolecuilderJobs_la_include_HEADERS = $(JOBSHEADER)
-libMolecuilderJobs_la_CPPFLAGS = $(AM_CPPFLAGS) ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} ${JobMarket_CFLAGS} -Dvmg_float=double -Dvmg_int=int $(VMG_CFLAGS)
-libMolecuilderJobs_Work_la_CPPFLAGS = $(AM_CPPFLAGS) ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} ${JobMarket_CFLAGS} -Dvmg_float=double -Dvmg_int=int $(VMG_CFLAGS)
+libMolecuilderJobs_la_CPPFLAGS = \
+	$(AM_CPPFLAGS) ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} ${JobMarket_CFLAGS} -Dvmg_float=double -Dvmg_int=int $(VMG_CFLAGS)
 libMolecuilderJobs_la_LDFLAGS = $(AM_LDFLAGS) \
   $(BOOST_IOSTREAMS_LDFLAGS) \
   $(BOOST_SERIALIZATION_LDFLAGS) \
   $(BOOST_SYSTEM_LDFLAGS)
-libMolecuilderJobs_Work_la_LDFLAGS = $(AM_LDFLAGS)
-if CONDJOBMARKET
 libMolecuilderJobs_la_LDFLAGS += \
   $(JobMarket_LDFLAGS)
-libMolecuilderJobs_Work_la_LDFLAGS += \
-  $(JobMarket_LDFLAGS)
-endif
-
 libMolecuilderJobs_la_LDFLAGS += \
   $(CodePatterns_LDFLAGS)
@@ -87,9 +119,10 @@
 endif
 endif
-
 if CONDJOBMARKET
 libMolecuilderJobs_la_LIBADD += \
   $(JobMarket_LIBS)
 libMolecuilderJobs_Work_la_LIBADD = \
+	$(top_builddir)/ThirdParty/mpqc_open/src/bin/mpqc/libmolecuilder_mpqc.la \
+	libmolecuilder_mpqc_extract.la \
   $(JobMarket_LIBS)
 endif
Index: src/Jobs/mpqc_extract.cc
===================================================================
--- src/Jobs/mpqc_extract.cc	(revision fbf005fe6f3de1a0fe7c1467b64db03e679102ae)
+++ src/Jobs/mpqc_extract.cc	(revision fbf005fe6f3de1a0fe7c1467b64db03e679102ae)
@@ -0,0 +1,606 @@
+//
+// mpqc_extract.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
+
+#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 "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
+
+#include <chemistry/qc/scf/linkage.h>
+#include <chemistry/qc/dft/linkage.h>
+#include <chemistry/qc/mbpt/linkage.h>
+#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_MBPTR12
+#  include <chemistry/qc/mbptr12/linkage.h>
+#endif
+#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_CINTS
+#  include <chemistry/qc/cints/linkage.h>
+#endif
+//#include <chemistry/qc/psi/linkage.h>
+#include <util/state/linkage.h>
+#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_CC
+#  include <chemistry/qc/cc/linkage.h>
+#endif
+#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_PSI
+#  include <chemistry/qc/psi/linkage.h>
+#endif
+#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_INTCCA
+#  include <chemistry/qc/intcca/linkage.h>
+#endif
+
+#include "mpqc_extract.h"
+
+using namespace std;
+using namespace sc;
+
+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;
+}
+
+/** 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)
+{
+  int 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;
+}
+
+/** Extractor function that is called after all calculations have been made.
+ *
+ * \param data result structure to fill
+ */
+void extractResults(
+	 Ref<MolecularEnergy> &mole,
+	 void *_data
+	 )
+{
+  MPQCData &data = *static_cast<MPQCData *>(_data);
+	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;
+}
+
+void extractTimings(
+	 Ref<RegionTimer> &tim,
+	 void *_data)
+{
+  MPQCData &data = *static_cast<MPQCData *>(_data);
+  // 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;
+}
+
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision 69682af4a24e4cabdca10e5bae32c1b59a630b37)
+++ src/Makefile.am	(revision fbf005fe6f3de1a0fe7c1467b64db03e679102ae)
@@ -429,5 +429,10 @@
 molecuilder_poolworker_SOURCES = poolworker.cpp 
 molecuilder_poolworker_LDFLAGS = $(AM_LDFLAGS) $(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LDFLAGS)
-molecuilder_poolworker_CPPFLAGS = $(AM_CPPFLAGS) $(JobMarket_CFLAGS)
+molecuilder_poolworker_CPPFLAGS = \
+	$(AM_CPPFLAGS) $(JobMarket_CFLAGS) \
+	-I$(top_builddir)/ThirdParty/mpqc_open/src/lib \
+	-I$(top_srcdir)/ThirdParty/mpqc_open/src/lib \
+	-I$(top_srcdir)/ThirdParty/mpqc_open/src/bin/mpqc
+
 molecuilder_poolworker_LDADD = \
 	libMolecuilderJobs.la \
Index: src/poolworker.cpp
===================================================================
--- src/poolworker.cpp	(revision 69682af4a24e4cabdca10e5bae32c1b59a630b37)
+++ src/poolworker.cpp	(revision fbf005fe6f3de1a0fe7c1467b64db03e679102ae)
@@ -36,14 +36,30 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include <string>
+
+#include <boost/bind.hpp>
+#include <boost/tokenizer.hpp>
+
 #include "JobMarket/poolworker_main.hpp"
+
+#include "mpqc.h"
 
 // force binding to libMolecuilderJobs
 // note we bind here to sort-of mock objects which are bases on SystemCommandJob
 #include "Jobs/MPQCCommandJob_binding.hpp"
-// poolworker has nothing to do with VMG jobs
+#include "Jobs/MPQCJob_binding.hpp"
+#include "Jobs/VMGJob_binding.hpp"
 
 int main(int argc, char* argv[])
 {
-  return poolworker_main(argc,argv);
+  mpqc::InitValues initvalues;
+
+  mpqc::init(initvalues, argc, argv);
+
+  const int exitflag = poolworker_main(argc, argv);
+
+  mpqc::cleanup(initvalues);
+
+  return exitflag;
 }
 
