source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 849cd8

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. 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/>.
 */

/*
 * FragmentationAutomationAction.cpp
 *
 *  Created on: May 18, 2012
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <boost/archive/text_iarchive.hpp>
// boost asio needs specific operator new
#include <boost/asio.hpp>

#include "CodePatterns/MemDebug.hpp"

#include <boost/mpl/remove.hpp>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "JobMarket/Controller/FragmentController.hpp"
#include "JobMarket/Jobs/FragmentJob.hpp"

#include "Atom/atom.hpp"
#include "Box.hpp"
#include "Fragmentation/EnergyMatrix.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Fragmentation/Fragmentation.hpp"
#include "Fragmentation/SetValues/Histogram.hpp"
#include "Fragmentation/SetValues/IndexedVectors.hpp"
#include "Fragmentation/HydrogenSaturation_enum.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
#include "Fragmentation/Summation/AllLevelOrthogonalSummator.hpp"
#include "Fragmentation/Summation/AllLevelSummator.hpp"
#include "Fragmentation/Summation/writeTable.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "Helpers/defs.hpp"
#include "Jobs/MPQCJob.hpp"
#include "Jobs/MPQCData.hpp"
#include "Jobs/MPQCData_printKeyNames.hpp"
#include "Jobs/Grid/SamplingGrid.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#ifdef HAVE_VMG
#include "Jobs/VMGJob.hpp"
#endif
#include "molecule.hpp"
#include "World.hpp"

#include <iostream>
#include <string>
#include <vector>

#include <boost/mpl/for_each.hpp>

#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"

using namespace MoleCuilder;

// and construct the stuff
#include "FragmentationAutomationAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */

class controller_AddOn;

// needs to be defined for using the FragmentController
controller_AddOn *getAddOn()
{
  return NULL;
}

const int LEVEL = 5;

/** Creates a MPQCCommandJob with argument \a filename.
 *
 * @param jobs created job is added to this vector
 * @param command mpqc command to execute
 * @param filename filename being argument to job
 * @param nextid id for this job
 */
void parsejob(
    std::vector<FragmentJob::ptr> &jobs,
    const std::string &command,
    const std::string &filename,
    const JobId_t nextid)
{
  std::ifstream file;
  file.open(filename.c_str());
  ASSERT( file.good(), "parsejob() - file "+filename+" does not exist.");
  std::string output((std::istreambuf_iterator<char>(file)),
      std::istreambuf_iterator<char>());
  double begin[NDIM] = { 0., 0., 0. };
  const RealSpaceMatrix& M = World::getInstance().getDomain().getM();
  const double size = M.at(0,0);
  ASSERT( M.determinant() == size*size*size,
      "parsejob() - current domain matrix "+toString(M)+" is not cubic.");
  const int level = LEVEL;
  FragmentJob::ptr testJob( new MPQCJob(nextid, output, begin, size, level) );
  jobs.push_back(testJob);
  file.close();
  LOG(1, "INFO: Added MPQCCommandJob from file "+filename+".");
}

/** Helper function to get number of atoms somehow.
 *
 * Here, we just parse the number of lines in the adjacency file as
 * it should correspond to the number of atoms, except when some atoms
 * are not bonded, but then fragmentation makes no sense.
 *
 * @param path path to the adjacency file
 */
size_t getNoAtomsFromAdjacencyFile(const std::string &path)
{
  size_t NoAtoms = 0;

  // parse in special file to get atom count (from line count)
  std::string filename(path);
  filename += FRAGMENTPREFIX;
  filename += ADJACENCYFILE;
  std::ifstream adjacency(filename.c_str());
  if (adjacency.fail()) {
    LOG(0, endl << "getNoAtomsFromAdjacencyFile() - Unable to open " << filename << ", is the directory correct?");
    return false;
  }
  std::string buffer;
  while (getline(adjacency, buffer))
    NoAtoms++;
  LOG(1, "INFO: There are " << NoAtoms << " atoms.");

  return NoAtoms;
}

/** Extracts MPQCData from received vector of FragmentResults.
 *
 * @param results results to extract MPQCData from
 * @param fragmentData on return array filled with extracted MPQCData
 */
void ConvertFragmentResultToMPQCData(
    const std::vector<FragmentResult::ptr> &results,
    std::vector<MPQCData> &fragmentData)
{
  // extract results
  fragmentData.clear();
  fragmentData.reserve(results.size());

  LOG(2, "DEBUG: Parsing now through " << results.size() << " results.");
  for (std::vector<FragmentResult::ptr>::const_iterator iter = results.begin();
      iter != results.end(); ++iter) {
    //LOG(1, "RESULT: job #"+toString((*iter)->getId())+": "+toString((*iter)->result));
    MPQCData extractedData;
    std::stringstream inputstream((*iter)->result);
    LOG(2, "DEBUG: First 50 characters FragmentResult's string: "+(*iter)->result.substr(0, 50));
    boost::archive::text_iarchive ia(inputstream);
    ia >> extractedData;
    LOG(1, "INFO: extracted data is " << extractedData << ".");
    fragmentData.push_back(extractedData);
  }

  ASSERT( results.size() == fragmentData.size(),
      "ConvertFragmentResultToMPQCData() - the number of extracted data differs from the number of results.");
}

/** Creates a lookup from FragmentJob::id to the true fragment  number.
 *
 * @param results result with job ids
 * @param MatrixNrLookup Lookup up-map, filled on return
 * @param FragmentCounter total number of fragments on return
 */
void createMatrixNrLookup(
    const std::vector<FragmentResult::ptr> &results,
    std::map< JobId_t, size_t > &MatrixNrLookup,
    size_t &FragmentCounter)
{
  // align fragments
  MatrixNrLookup.clear();
  FragmentCounter = 0;
  {
    // bring ids in order ...
    typedef std::map< JobId_t, FragmentResult::ptr> IdResultMap_t;
    IdResultMap_t IdResultMap;
    for (std::vector<FragmentResult::ptr>::const_iterator iter = results.begin();
        iter != results.end(); ++iter) {
  #ifndef NDEBUG
      std::pair< IdResultMap_t::iterator, bool> inserter =
  #endif
      IdResultMap.insert( make_pair((*iter)->getId(), *iter) );
      ASSERT( inserter.second,
          "ExtractMPQCDataFromResults() - two results have same id "
          +toString((*iter)->getId())+".");
    }
    // ... and fill lookup
    for(IdResultMap_t::const_iterator iter = IdResultMap.begin();
        iter != IdResultMap.end(); ++iter)
      MatrixNrLookup.insert( make_pair(iter->first, FragmentCounter++) );
  }
  LOG(1, "INFO: There are " << FragmentCounter << " fragments.");
}

/** Place results from FragmentResult into EnergyMatrix and ForceMatrix.
 *
 * @param results results with ids to associate with fragment number
 * @param fragmentData MPQCData resulting from the jobs
 * @param MatrixNrLookup Lookup up-map from job id to fragment number
 * @param FragmentCounter total number of fragments
 * @param NoAtoms total number of atoms
 * @param Energy energy matrix to be filled on return
 * @param Force force matrix to be filled on return
 * @return true - everything ok, false - else
 */
bool putResultsintoMatrices(
    const std::vector<FragmentResult::ptr> &results,
    const std::vector<MPQCData> &fragmentData,
    std::map< JobId_t, size_t > &MatrixNrLookup,
    const size_t FragmentCounter,
    const size_t NoAtoms,
    EnergyMatrix &Energy,
    ForceMatrix &Force)
{
  ASSERT( results.size() == fragmentData.size(),
      "printReceivedMPQCResults() - results and fragmentData differ in size.");
  std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
  std::vector<FragmentResult::ptr>::const_iterator resultiter = results.begin();
  for (; dataiter != fragmentData.end(); ++dataiter, ++resultiter) {
    const MPQCData &extractedData = *dataiter;
    // place results into EnergyMatrix ...
    {
      MatrixContainer::MatrixArray matrix;
      matrix.resize(1);
      matrix[0].resize(1, extractedData.energies.total);
      if (!Energy.AddMatrix(
          std::string("MPQCJob ")+toString((*resultiter)->getId()),
          matrix,
          MatrixNrLookup[(*resultiter)->getId()])) {
        ELOG(1, "Adding energy matrix failed.");
        return false;
      }
    }
    // ... and ForceMatrix (with two empty columns in front)
    {
      MatrixContainer::MatrixArray matrix;
      const size_t rows = extractedData.forces.size();
      matrix.resize(rows);
      for (size_t i=0;i<rows;++i) {
        const size_t columns = 2+extractedData.forces[i].size();
        matrix[i].resize(columns, 0.);
  //      for (size_t j=0;j<2;++j)
  //        matrix[i][j] = 0.;
        for (size_t j=2;j<columns;++j)
          matrix[i][j] = extractedData.forces[i][j-2];
      }
      if (!Force.AddMatrix(
          std::string("MPQCJob ")+toString((*resultiter)->getId()),
          matrix,
          MatrixNrLookup[(*resultiter)->getId()])) {
        ELOG(1, "Adding force matrix failed.");
        return false;
      }
    }
  }
  // add one more matrix (not required for energy)
  MatrixContainer::MatrixArray matrix;
  matrix.resize(1);
  matrix[0].resize(1, 0.);
  if (!Energy.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
    return false;
  // but for energy because we need to know total number of atoms
  matrix.resize(NoAtoms);
  for (size_t i = 0; i< NoAtoms; ++i)
    matrix[i].resize(2+NDIM, 0.);
  if (!Force.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
    return false;

  return true;
}

template <typename T>
void convertMPQCDataTo(
    const std::vector<MPQCData> &fragmentData,
    std::vector<T> &MPQCData_fused)
{
  MPQCData_fused.clear();
}

template <>
void convertMPQCDataTo<MPQCDataEnergyMap_t>(
    const std::vector<MPQCData> &fragmentData,
    std::vector<MPQCDataEnergyMap_t> &MPQCData_Energy_fused)
{
  // energy_t
  MPQCData_Energy_fused.clear();
  MPQCData_Energy_fused.reserve(fragmentData.size());
  for(std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
      dataiter != fragmentData.end(); ++dataiter) {
    const MPQCData &extractedData = *dataiter;
    LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
    MPQCDataEnergyMap_t instance;
    boost::fusion::at_key<MPQCDataFused::energy_total>(instance) = extractedData.energies.total;
    boost::fusion::at_key<MPQCDataFused::energy_nuclear_repulsion>(instance) = extractedData.energies.nuclear_repulsion;
    boost::fusion::at_key<MPQCDataFused::energy_electron_coulomb>(instance) = extractedData.energies.electron_coulomb;
    boost::fusion::at_key<MPQCDataFused::energy_electron_exchange>(instance) = extractedData.energies.electron_exchange;
    boost::fusion::at_key<MPQCDataFused::energy_correlation>(instance) = extractedData.energies.correlation;
    boost::fusion::at_key<MPQCDataFused::energy_overlap>(instance) = extractedData.energies.overlap;
    boost::fusion::at_key<MPQCDataFused::energy_kinetic>(instance) = extractedData.energies.kinetic;
    boost::fusion::at_key<MPQCDataFused::energy_hcore>(instance) = extractedData.energies.hcore;
    boost::fusion::at_key<MPQCDataFused::energy_eigenvalues>(instance) = extractedData.energies.eigenvalues;
    MPQCData_Energy_fused.push_back(instance);
  }
}


void convertMPQCDatatoForceMap(
    const std::vector<MPQCData> &fragmentData,
    const KeySetsContainer &ForceKeySet,
    std::vector<MPQCDataForceMap_t> &MPQCData_Force_fused)
{
  // forces
  ASSERT( ForceKeySet.KeySets.size() == fragmentData.size(),
      "FragmentationAutomationAction::performCall() - indices and fragmentData differ in size.");
  MPQCData_Force_fused.clear();
  MPQCData_Force_fused.reserve(fragmentData.size());
  std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
  KeySetsContainer::ArrayOfIntVectors::const_iterator arrayiter = ForceKeySet.KeySets.begin();
  for(;dataiter != fragmentData.end(); ++dataiter, ++arrayiter) {
    const MPQCData &extractedData = *dataiter;
    LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
    MPQCDataForceMap_t instance;
    // must convert int to index_t
    IndexedVectors::indices_t indices(arrayiter->begin(), arrayiter->end());
    boost::fusion::at_key<MPQCDataFused::forces>(instance) =
        IndexedVectors(indices, extractedData.forces);
    MPQCData_Force_fused.push_back(instance);
  }
}

template <>
void convertMPQCDataTo<MPQCDataGridMap_t>(
    const std::vector<MPQCData> &fragmentData,
    std::vector<MPQCDataGridMap_t> &MPQCData_Grid_fused)
{
  // sampled_grid
  MPQCData_Grid_fused.clear();
  MPQCData_Grid_fused.reserve(fragmentData.size());
  for(std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
      dataiter != fragmentData.end(); ++dataiter) {
    const MPQCData &extractedData = *dataiter;
    LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
    MPQCDataGridMap_t instance;
    boost::fusion::at_key<MPQCDataFused::sampled_grid>(instance) = extractedData.sampled_grid;
    MPQCData_Grid_fused.push_back(instance);
  }
}

template <>
void convertMPQCDataTo<MPQCDataTimeMap_t>(
    const std::vector<MPQCData> &fragmentData,
    std::vector<MPQCDataTimeMap_t> &MPQCData_Time_fused)
{
  // times
  MPQCData_Time_fused.clear();
  MPQCData_Time_fused.reserve(fragmentData.size());
  for(std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
      dataiter != fragmentData.end(); ++dataiter) {
    const MPQCData &extractedData = *dataiter;
    LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
    MPQCDataTimeMap_t instance;
    boost::fusion::at_key<MPQCDataFused::times_walltime>(instance) = extractedData.times.walltime;
    boost::fusion::at_key<MPQCDataFused::times_cputime>(instance) = extractedData.times.cputime;
    boost::fusion::at_key<MPQCDataFused::times_flops>(instance) = extractedData.times.flops;
    MPQCData_Time_fused.push_back(instance);
  }
}

template <typename TypeMap, typename TypeVector>
std::vector<TypeMap> SumUpPerLevel(
    const std::vector<MPQCData> &fragmentData,
    const std::vector<JobId_t> &jobids,
    std::map< JobId_t, size_t > &MatrixNrLookup,
    const IndexSetContainer::ptr &container,
    SubsetMap::ptr &subsetmap
    )
{
  // place data into boost::fusion::map instance
  std::vector<TypeMap> MPQCData_fused;
  convertMPQCDataTo<TypeMap>(fragmentData, MPQCData_fused);
  // instantiate summator
  std::vector<TypeMap> Result_fused(subsetmap->getMaximumSubsetLevel());
  AllLevelOrthogonalSummator<TypeMap> Summer(
              subsetmap,
              MPQCData_fused,
              jobids,
              container->getContainer(),
              MatrixNrLookup,
              Result_fused);
  // sum up for each type key in TypeVector
  boost::mpl::for_each<TypeVector>(boost::ref(Summer));
  return Result_fused;
}

/** Print MPQCData from received results.
 *
 * @param results results with ids to associate with fragment number
 * @param fragmentData MPQCData resulting from the jobs
 * @param KeySetFilename filename with keysets to associate forces correctly
 * @param NoAtoms total number of atoms
 * @param full_sample summed up charge from fragments on return
 */
bool sumUpChargeDensity(
    const std::vector<FragmentResult::ptr> &results,
    const std::vector<MPQCData> &fragmentData,
    const std::string &KeySetFilename,
    SamplingGrid &full_sample)
{
  // create lookup from job nr to fragment number
  std::map< JobId_t, size_t > MatrixNrLookup;
  size_t FragmentCounter = 0;
  createMatrixNrLookup(results, MatrixNrLookup, FragmentCounter);

  // initialise keysets
  KeySetsContainer KeySet;
  {
    // else needs keysets without hydrogens
    std::stringstream filename;
    filename << FRAGMENTPREFIX << KEYSETFILE;
    if (!KeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
  }

  /// prepare for OrthogonalSummation

  // convert KeySetContainer to IndexSetContainer
  IndexSetContainer::ptr container(new IndexSetContainer(KeySet));
  // create the map of all keysets
  SubsetMap::ptr subsetmap(new SubsetMap(*container));

  // create a vector of all job ids
  std::vector<JobId_t> jobids(results.size(), JobId::IllegalJob);
  std::transform(results.begin(), results.end(), jobids.begin(),
      boost::bind(&FragmentResult::getId,
          boost::bind(&FragmentResult::ptr::operator->, _1)));

  /// convert all MPQCData to MPQCDataMap_t
  std::vector<MPQCDataGridMap_t> Result_Grid_fused(
    SumUpPerLevel<MPQCDataGridMap_t, MPQCDataGridVector_t>(
        fragmentData, jobids, MatrixNrLookup, container, subsetmap));
  // obtain full grid
  full_sample = boost::fusion::at_key<MPQCDataFused::sampled_grid>(Result_Grid_fused.back());

  return true;
}


/** Print MPQCData from received results.
 *
 * @param results results with ids to associate with fragment number
 * @param fragmentData MPQCData resulting from the jobs
 * @param KeySetFilename filename with keysets to associate forces correctly
 * @param NoAtoms total number of atoms
 * @param full_sample summed up charge from fragments on return
 */
bool printReceivedMPQCResults(
    const std::vector<FragmentResult::ptr> &results,
    const std::vector<MPQCData> &fragmentData,
    const std::string &KeySetFilename,
    size_t NoAtoms,
    SamplingGrid &full_sample)
{
  // create lookup from job nr to fragment number
  std::map< JobId_t, size_t > MatrixNrLookup;
  size_t FragmentCounter = 0;
  createMatrixNrLookup(results, MatrixNrLookup, FragmentCounter);

  // place results into maps
  EnergyMatrix Energy;
  ForceMatrix Force;
  if (!putResultsintoMatrices(results, fragmentData, MatrixNrLookup, FragmentCounter, NoAtoms, Energy, Force))
    return false;

  // initialise keysets
  KeySetsContainer KeySet;
  KeySetsContainer ForceKeySet;
  if (!Energy.InitialiseIndices()) return false;

  if (!Force.ParseIndices(KeySetFilename.c_str())) return false;

  {
    // else needs keysets without hydrogens
    std::stringstream filename;
    filename << FRAGMENTPREFIX << KEYSETFILE;
    if (!KeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
  }

  {
    // forces need keysets including hydrogens
    std::stringstream filename;
    filename << FRAGMENTPREFIX << FORCESFILE;
    if (!ForceKeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
  }

  /// prepare for OrthogonalSummation

  // convert KeySetContainer to IndexSetContainer
  IndexSetContainer::ptr container(new IndexSetContainer(KeySet));
  // create the map of all keysets
  SubsetMap::ptr subsetmap(new SubsetMap(*container));

  // create a vector of all job ids
  std::vector<JobId_t> jobids(results.size(), JobId::IllegalJob);
  std::transform(results.begin(), results.end(), jobids.begin(),
      boost::bind(&FragmentResult::getId,
          boost::bind(&FragmentResult::ptr::operator->, _1)));

  /// convert all MPQCData to MPQCDataMap_t
  {
    ASSERT( ForceKeySet.KeySets.size() == fragmentData.size(),
        "FragmentationAutomationAction::performCall() - ForceKeySet's KeySets and fragmentData differ in size.");

    std::vector<MPQCDataEnergyMap_t> Result_Energy_fused(
        SumUpPerLevel<MPQCDataEnergyMap_t, MPQCDataEnergyVector_t>(
            fragmentData, jobids, MatrixNrLookup, container, subsetmap));
    std::vector<MPQCDataGridMap_t> Result_Grid_fused(
      SumUpPerLevel<MPQCDataGridMap_t, MPQCDataGridVector_t>(
          fragmentData, jobids, MatrixNrLookup, container, subsetmap));
    std::vector<MPQCDataTimeMap_t> Result_Time_fused(
      SumUpPerLevel<MPQCDataTimeMap_t, MPQCDataTimeVector_t>(
          fragmentData, jobids, MatrixNrLookup, container, subsetmap));

    // force has extra converter
    std::vector<MPQCDataForceMap_t> MPQCData_Force_fused;
    convertMPQCDatatoForceMap(fragmentData, ForceKeySet, MPQCData_Force_fused);
    std::vector<MPQCDataForceMap_t> Result_Force_fused(subsetmap->getMaximumSubsetLevel());
    AllLevelOrthogonalSummator<MPQCDataForceMap_t> forceSummer(
                subsetmap,
                MPQCData_Force_fused,
                jobids,
                container->getContainer(),
                MatrixNrLookup,
                Result_Force_fused);
    boost::mpl::for_each<MPQCDataForceVector_t>(boost::ref(forceSummer));

    // obtain full grid
    full_sample = boost::fusion::at_key<MPQCDataFused::sampled_grid>(Result_Grid_fused.back());

    // print tables (without eigenvalues, they go extra)
    const size_t MaxLevel = subsetmap->getMaximumSubsetLevel();
    typedef boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type MPQCDataEnergyVector_noeigenvalues_t;
    const std::string energyresult =
        writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
            Result_Energy_fused, MaxLevel);
    LOG(0, "Energy table is \n" << energyresult);
    const std::string eigenvalueresult;

    LOG(0, "Eigenvalue table is \n" << eigenvalueresult);
    const std::string forceresult =
        writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
            Result_Force_fused, MaxLevel);
    LOG(0, "Force table is \n" << forceresult);
    // we don't want to print grid to a table
    // print times (without flops for now)
    typedef boost::mpl::remove<MPQCDataTimeVector_t, MPQCDataFused::times_flops>::type MPQCDataTimeVector_noflops_t;
    const std::string timesresult =
        writeTable<MPQCDataTimeMap_t, MPQCDataTimeVector_noflops_t >()(
            Result_Time_fused, MaxLevel);
    LOG(0, "Times table is \n" << timesresult);
  }

  // combine all found data
  if (!KeySet.ParseManyBodyTerms()) return false;

  EnergyMatrix EnergyFragments;
  ForceMatrix ForceFragments;
  if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return false;
  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return false;

  if(!Energy.SetLastMatrix(0., 0)) return false;
  if(!Force.SetLastMatrix(0., 2)) return false;

  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    // --------- sum up energy --------------------
    LOG(1, "INFO: Summing energy of order " << BondOrder+1 << " ...");
    if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return false;
    if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return false;

    // --------- sum up Forces --------------------
    LOG(1, "INFO: Summing forces of order " << BondOrder+1 << " ...");
    if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return false;
    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return false;
  }

  // for debugging print resulting energy and forces
  LOG(1, "INFO: Resulting energy is " << Energy.Matrix[ FragmentCounter ][0][0]);
  std::stringstream output;
  for (int i=0; i< Force.RowCounter[FragmentCounter]; ++i) {
    for (int j=0; j< Force.ColumnCounter[FragmentCounter]; ++j)
      output << Force.Matrix[ FragmentCounter ][i][j] << " ";
    output << "\n";
  }
  LOG(1, "INFO: Resulting forces are " << std::endl << output.str());

  return true;
}


void RunService(
    boost::asio::io_service &io_service,
    std::string message)
{
  message = std::string("io_service: ") + message;
  io_service.reset();
  Info info(message.c_str());
  io_service.run();
}

void requestIds(
    FragmentController &controller,
    const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
    const size_t numberjobs)
{
  controller.requestIds(params.host.get(), params.port.get(), numberjobs);
}

bool createJobsFromFiles(
    FragmentController &controller,
    const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
    const std::vector< boost::filesystem::path > &jobfiles)
{
  std::vector<FragmentJob::ptr> jobs;
  for (std::vector< boost::filesystem::path >::const_iterator iter = jobfiles.begin();
      iter != jobfiles .end(); ++iter) {
    const std::string &filename = (*iter).string();
    if (boost::filesystem::exists(filename)) {
      const JobId_t next_id = controller.getAvailableId();
      LOG(1, "INFO: Creating MPQCCommandJob with filename'"
          +filename+"', and id "+toString(next_id)+".");
      parsejob(jobs, params.executable.get().string(), filename, next_id);
    } else {
      ELOG(1, "Fragment job "+filename+" does not exist.");
      return false;
    }
  }
  controller.addJobs(jobs);
  controller.sendJobs(params.host.get(), params.port.get());
  return true;
}

#ifdef HAVE_VMG
bool createLongRangeJobs(
    FragmentController &controller,
    const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
    const std::vector<MPQCData> &fragmentData,
    const SamplingGrid &full_sampled_grid)
{
  std::vector<FragmentJob::ptr> jobs;
  // add one job for each fragment as the short-range correction which we need
  // to subtract from the obtained full potential to get the long-range part only
  for (std::vector<MPQCData>::const_iterator iter = fragmentData.begin();
      iter != fragmentData.end(); ++iter) {
    const JobId_t next_id = controller.getAvailableId();
    LOG(1, "INFO: Creating VMGJob with " << iter->sampled_grid.sampled_grid.size()
        << " gridpoints and " << iter->charges.size() << " particle charges.");
    FragmentJob::ptr testJob(
        new VMGJob(next_id, iter->sampled_grid, iter->positions, iter->charges) );
    jobs.push_back(testJob);
  }

  {
//    const int level = LEVEL;
//    const int GRID = pow(2, level);
//    double begin[NDIM] = { 0., 0., 0. };
//    const RealSpaceMatrix& M = World::getInstance().getDomain().getM();
//    const double size = M.at(0,0);
//    ASSERT( M.determinant() == size*size*size,
//        "createLongRangeJobs() - current domain matrix "+toString(M)+" is not cubic.");
    const std::vector< std::vector<double> > positions;
    const std::vector<double> charges;
    const JobId_t next_id = controller.getAvailableId();
    LOG(1, "INFO: Creating full VMGJob with " << full_sampled_grid.sampled_grid.size()
        << " gridpoints and " << charges.size() << " particle charges.");
    FragmentJob::ptr testJob(
        new VMGJob(next_id, full_sampled_grid, positions, charges) );
    jobs.push_back(testJob);
  }

  // then send jobs to controller
  controller.addJobs(jobs);
  controller.sendJobs(params.host.get(), params.port.get());
  return true;
}
#endif

void WaitforResults(
    boost::asio::io_service &io_service,
    FragmentController &controller,
    const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
    const size_t NoExpectedResults
    )
{
  size_t NoCalculatedResults = 0;
  while (NoCalculatedResults != NoExpectedResults) {
    // wait a bit
    boost::asio::deadline_timer timer(io_service);
    timer.expires_from_now(boost::posix_time::milliseconds(500));
    timer.wait();
    // then request status
    controller.checkResults(params.host.get(), params.port.get());
    RunService(io_service, "Checking on results");

    const std::pair<size_t, size_t> JobStatus = controller.getJobStatus();
    LOG(1, "INFO: #" << JobStatus.first << " are waiting in the queue and #" << JobStatus.second << " jobs are calculated so far.");
    NoCalculatedResults = JobStatus.second;
  }
}


Action::state_ptr FragmentationFragmentationAutomationAction::performCall() {
  boost::asio::io_service io_service;
  FragmentController controller(io_service);

  // TODO: Have io_service run in second thread and merge with current again eventually

  // Phase One: obtain ids
  std::vector< boost::filesystem::path > jobfiles = params.jobfiles.get();
  requestIds(controller, params, jobfiles.size());
  RunService(io_service, "Requesting ids");

  // Phase Two: create and add MPQCJobs
  if (!createJobsFromFiles(controller, params, jobfiles))
    return Action::failure;
  RunService(io_service, "Adding MPQCJobs");

  // Phase Three: calculate result
  WaitforResults(io_service, controller, params, jobfiles.size());
  controller.receiveResults(params.host.get(), params.port.get());
  RunService(io_service, "Requesting short-range results");
  std::vector<FragmentResult::ptr> MPQCresults = controller.getReceivedResults();
  std::vector<MPQCData> fragmentData;
  ConvertFragmentResultToMPQCData(MPQCresults, fragmentData);

#ifdef HAVE_VMG
  if (params.DoLongrange.get()) {
  // obtain combined charge density
  LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
  SamplingGrid full_sample;
  sumUpChargeDensity(
      MPQCresults,
      fragmentData,
      params.path.get(),
      full_sample);

  // Phase Four: obtain more ids
  requestIds(controller, params, fragmentData.size()+1);
  RunService(io_service, "Requesting ids");

  // Phase Five: create VMGJobs
  if (!createLongRangeJobs(controller, params, fragmentData, full_sample))
    return Action::failure;
  RunService(io_service, "Adding VMGJobs");

  // Phase Six: calculate result
  WaitforResults(io_service, controller, params, fragmentData.size()+1);
  controller.receiveResults(params.host.get(), params.port.get());
  RunService(io_service, "Requesting long-range results");
  std::vector<FragmentResult::ptr> VMGresults = controller.getReceivedResults();
  ASSERT( MPQCresults.size()+1 == VMGresults.size(),
      "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresultd and VMGresults don't match.");

  // Final phase: print result
  {
    LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
    printReceivedMPQCResults(
        MPQCresults,
        fragmentData,
        params.path.get(),
        getNoAtomsFromAdjacencyFile(params.path.get()),
        full_sample);
  }
  }
#endif
  size_t Exitflag = controller.getExitflag();

  return (Exitflag == 0) ? Action::success : Action::failure;
}

Action::state_ptr FragmentationFragmentationAutomationAction::performUndo(Action::state_ptr _state) {
  return Action::success;
}

Action::state_ptr FragmentationFragmentationAutomationAction::performRedo(Action::state_ptr _state){
  return Action::success;
}

bool FragmentationFragmentationAutomationAction::canUndo() {
  return false;
}

bool FragmentationFragmentationAutomationAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */