source: src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp@ eaef9b

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

/*
 * AnalyseFragmentationResultsAction.cpp
 *
 *  Created on: Mar 8, 2013
 *      Author: heber
 */

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

// include headers that implement a archive in simple text format
// and before MemDebug due to placement new
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include "CodePatterns/MemDebug.hpp"

#include <boost/foreach.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/mpl/remove.hpp>

#include <algorithm>
#include <fstream>
#include <iostream>
//#include <numeric>
#include <string>
#include <vector>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

#ifdef HAVE_JOBMARKET
#include "JobMarket/types.hpp"
#else
typedef size_t JobId_t;
#endif

#include "Descriptors/AtomIdDescriptor.hpp"
#include "Atom/atom.hpp"
#include "Element/element.hpp"
#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
#include "Fragmentation/Summation/Containers/FragmentationChargeDensity.hpp"
#include "Fragmentation/Summation/Containers/FragmentationResultContainer.hpp"
#include "Fragmentation/Summation/Containers/FragmentationShortRangeResults.hpp"
#include "Fragmentation/Summation/Containers/MPQCData.hpp"
#include "Fragmentation/Summation/Containers/MPQCData_printKeyNames.hpp"
#include "Fragmentation/Homology/HomologyContainer.hpp"
#include "Fragmentation/Homology/HomologyGraph.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
#include "Fragmentation/Summation/SetValues/Eigenvalues.hpp"
#include "Fragmentation/Summation/SetValues/Fragment.hpp"
#include "Fragmentation/Summation/SetValues/FragmentForces.hpp"
#include "Fragmentation/Summation/SetValues/Histogram.hpp"
#include "Fragmentation/Summation/SetValues/IndexedVectors.hpp"
#include "Fragmentation/Summation/IndexSetContainer.hpp"
#include "Fragmentation/Summation/writeIndexedTable.hpp"
#include "Fragmentation/Summation/writeTable.hpp"
#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
#include "Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp"
#include "Fragmentation/Summation/Containers/GridDownsampler.hpp"
#include "Fragmentation/Summation/Containers/VMGData.hpp"
#include "Fragmentation/Summation/Containers/VMGDataFused.hpp"
#include "Fragmentation/Summation/Containers/VMGDataMap.hpp"
#include "Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp"
#endif
#include "Helpers/defs.hpp"
#include "Potentials/Particles/ParticleRegistry.hpp"
#include "World.hpp"

#include "Actions/FragmentationAction/AnalyseFragmentationResultsAction.hpp"

using namespace MoleCuilder;

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

void writeToFile(const std::string &filename, const std::string contents)
{
  std::ofstream tablefile(filename.c_str());
  tablefile << contents;
  tablefile.close();
}

/** Print cycle correction from received results.
 *
 * @param results summed up results container
 */
void printReceivedCycleResults(
    const FragmentationShortRangeResults &results)
{
  typedef boost::mpl::remove<
        boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type,
      MPQCDataFused::energy_eigenhistogram>::type
    MPQCDataEnergyVector_noeigenvalues_t;
  const std::string energyresult =
      writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
          results.Result_Energy_fused, results.getMaxLevel());
  LOG(2, "DEBUG: Energy table is \n" << energyresult);
  std::string filename;
  filename += FRAGMENTPREFIX + std::string("_CycleEnergy.dat");
  writeToFile(filename, energyresult);
}

/** Print (short range) energy, forces, and timings from received results per index set.
 *
 * @param results summed up results container
 */
void printReceivedShortResultsPerIndex(
    const FragmentationShortRangeResults &results)
{
  // print tables  per keyset(without eigenvalues, they go extra)
  typedef boost::mpl::remove<
        boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type,
      MPQCDataFused::energy_eigenhistogram>::type
    MPQCDataEnergyVector_noeigenvalues_t;
  const std::string energyresult =
      writeIndexedTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
          results.Result_perIndexSet_Energy, results.getMaxLevel());
  LOG(2, "DEBUG: Indexed Energy table is \n" << energyresult);
  std::string filename;
  filename += FRAGMENTPREFIX + std::string("_IndexedEnergy.dat");
  writeToFile(filename, energyresult);
}

/** Print (short range) energy, forces, and timings from received results.
 *
 * @param results summed up results container
 */
void printReceivedShortResults(
    const FragmentationShortRangeResults &results)
{
  // print tables (without eigenvalues, they go extra)
  {
    typedef boost::mpl::remove<
          boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type,
        MPQCDataFused::energy_eigenhistogram>::type
      MPQCDataEnergyVector_noeigenvalues_t;
    const std::string energyresult =
        writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
            results.Result_Energy_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Energy table is \n" << energyresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Energy.dat");
    writeToFile(filename, energyresult);
  }

  {
    typedef boost::mpl::list<
        MPQCDataFused::energy_eigenvalues
    > MPQCDataEigenvalues_t;
    const std::string eigenvalueresult =
        writeTable<MPQCDataEnergyMap_t, MPQCDataEigenvalues_t >()(
            results.Result_Energy_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Eigenvalue table is \n" << eigenvalueresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Eigenvalues.dat");
    writeToFile(filename, eigenvalueresult);
  }

  {
    typedef boost::mpl::list<
        MPQCDataFused::energy_eigenhistogram
    > MPQCDataEigenhistogram_t;
    const std::string eigenhistogramresult =
        writeTable<MPQCDataEnergyMap_t, MPQCDataEigenhistogram_t >()(
            results.Result_Energy_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Eigenhistogram table is \n" << eigenhistogramresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Eigenhistogram.dat");
    writeToFile(filename, eigenhistogramresult);
  }

  {
    typedef boost::mpl::list<
        MPQCDataFused::energy_eigenvalues
    > MPQCDataEigenvalues_t;
    const std::string eigenvalueresult =
        writeTable<MPQCDataEnergyMap_t, MPQCDataEigenvalues_t >()(
            results.Result_Energy_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Eigenvalue table is \n" << eigenvalueresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Eigenvalues.dat");
    writeToFile(filename, eigenvalueresult);
  }

  {
    const std::string forceresult =
        writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
            results.Result_Force_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Force table is \n" << forceresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Forces.dat");
    writeToFile(filename, forceresult);
  }
  // we don't want to print grid to a table
  {
    // print times (without flops for now)
    typedef boost::mpl::remove<
        boost::mpl::remove<MPQCDataTimeVector_t, MPQCDataFused::times_total_flops>::type,
        MPQCDataFused::times_gather_flops>::type
        MPQCDataTimeVector_noflops_t;
    const std::string timesresult =
        writeTable<MPQCDataTimeMap_t, MPQCDataTimeVector_noflops_t >()(
            results.Result_Time_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Times table is \n" << timesresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Times.dat");
    writeToFile(filename, timesresult);
  }
}


#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
/** Print long range energy from received results.
 *
 * @param results summed up results container
 */
void printReceivedFullResults(
    const FragmentationLongRangeResults &results)
{
  // print tables per keyset(without grids, they go extra)

  {
    const std::string gridresult =
        writeTable<VMGDataMap_t, VMGDataVector_t >()(
            results.Result_LongRange_fused, results.getMaxLevel(), 1);
    LOG(2, "DEBUG: VMG table is \n" << gridresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_VMGEnergy.dat");
    writeToFile(filename, gridresult);
  }

  if (results.hasLongRangeForces()) {
    const std::string forceresult =
        writeTable<VMGDataForceMap_t, VMGDataForceVector_t>()(
            results.Result_ForceLongRange_fused, results.getMaxLevel());
    LOG(2, "DEBUG: Force table is \n" << forceresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_VMGForces.dat");
    writeToFile(filename, forceresult);
  }

  {
    const std::string gridresult =
        writeTable<VMGDataLongRangeMap_t, VMGDataLongRangeVector_t >()(
            results.Result_LongRangeIntegrated_fused, results.getMaxLevel(), 1);
    LOG(2, "DEBUG: LongRange table is \n" << gridresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_LongRangeEnergy.dat");
    writeToFile(filename, gridresult);
  }

  if (results.hasLongRangeForces()) {
    const std::string forceresult =
        writeTable<VMGDataLongRangeForceMap_t, VMGDataLongRangeForceVector_t >()(
            results.Result_ForcesLongRangeIntegrated_fused, results.getMaxLevel(), 1);
    LOG(2, "DEBUG: ForcesLongRange table is \n" << forceresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_LongRangeForces.dat");
    writeToFile(filename, forceresult);
  }
}
#endif

bool appendToHomologies(
    const FragmentationShortRangeResults &shortrangeresults,
#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    const FragmentationLongRangeResults &longrangeresults,
#endif
    const bool storeGrids
    )
{
  /// read homology container (if present)
  HomologyContainer &homology_container = World::getInstance().getHomologies();

  /// append all fragments to a HomologyContainer
  HomologyContainer::container_t values;

  // convert KeySetContainer to IndexSetContainer
  // BUG: Conversion changes order w.r.t to Indices(!)
  const IndexSetContainer::Container_t &container = shortrangeresults.getContainer();
  const KeySetsContainer &Indices = shortrangeresults.getKeySet();
  const KeySetsContainer &ForceIndices = shortrangeresults.getForceKeySet();
  if (Indices.KeySets.size() != ForceIndices.KeySets.size()) {
    ELOG(1, "appendToHomologies() - Indices (" << Indices.KeySets.size()
        << ") and ForceIndices (" << ForceIndices.KeySets.size() << ") sizes differ.");
    return false;
  }
  KeySetsContainer::ArrayOfIntVectors::const_iterator iter = Indices.KeySets.begin();
  KeySetsContainer::ArrayOfIntVectors::const_iterator forceiter = ForceIndices.KeySets.begin();

  /// go through all fragments
  for (;iter != Indices.KeySets.end(); ++iter, ++forceiter) // go through each IndexSet
  {
    /// create new graph entry in HomologyContainer which is (key,value) type
    LOG(1, "INFO: Creating new graph with " << *forceiter << ".");
    IndexSet forceindexset;
    forceindexset.insert(forceiter->begin(), forceiter->end());
    const HomologyGraph graph(forceindexset);
    // get index set entity from container for lookup to fragment result
    IndexSet::ptr indexset(new IndexSet);
    indexset->insert(iter->begin(), iter->end());
    if (!forceindexset.contains(*indexset)) {
      // this caught an error with faulty KeySetsContainer::insert().
      // indexset and forceindexset need to be in same order and differ
      // only in forceindexset contains extra indices for saturation hydrogens
      ELOG(1, "appendToHomologies() - force index set " << forceindexset
          << " does not contain index set " << (*indexset) << ".");
      return false;
    }
    const IndexSetContainer::Container_t::const_iterator indexiter =
        std::lower_bound(container.begin(), container.end(), indexset);
    const IndexSet::ptr &index = *indexiter;
    ASSERT( *index == *indexset,
        "appendToHomologies() - could not find index set "+toString(*indexset)
        +" in the sorted IndexSetContainer.");
 

    /// we fill the value structure
    HomologyContainer::value_t value;
    value.containsGrids = storeGrids;
    // obtain fragment
    std::map<IndexSet::ptr, std::pair< MPQCDataFragmentMap_t, MPQCDataFragmentMap_t> >::const_iterator fragmentiter
        = shortrangeresults.Result_perIndexSet_Fragment.find(index);
    if (fragmentiter == shortrangeresults.Result_perIndexSet_Fragment.end()) {
      ELOG(1, "appendToHomologyFile() - cannot find index " << (*index)
          << " in FragmentResults.");
      return false;
    }
    value.fragment = boost::fusion::at_key<MPQCDataFused::fragment>(fragmentiter->second.first);
    LOG(1, "DEBUG: Created graph " << graph << " to fragment " << value.fragment << ".");

    // obtain energy
    std::map<IndexSet::ptr, std::pair< MPQCDataEnergyMap_t, MPQCDataEnergyMap_t> >::const_iterator energyiter
        = shortrangeresults.Result_perIndexSet_Energy.find(index);
    if (energyiter == shortrangeresults.Result_perIndexSet_Energy.end()) {
      ELOG(1, "appendToHomologyFile() - cannot find index " << (*index)
          << " in FragmentResults.");
      return false;
    }
    value.energy = boost::fusion::at_key<MPQCDataFused::energy_total>(energyiter->second.second); // contributions

    // only store sampled grids if desired
    if (storeGrids) {
#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
      // obtain charge distribution
      std::map<IndexSet::ptr, std::pair< MPQCDataGridMap_t, MPQCDataGridMap_t> >::const_iterator chargeiter
          = longrangeresults.Result_perIndexSet_Grid.find(index);
      if (chargeiter == longrangeresults.Result_perIndexSet_Grid.end()) {
        ELOG(1, "appendToHomologyFile() - cannot find index " << (*index)
            << " in FragmentResults.");
        return false;
      }
      value.charge_distribution = boost::fusion::at_key<MPQCDataFused::sampled_grid>(chargeiter->second.second); // contributions

      // obtain potential distribution
      std::map<IndexSet::ptr, std::pair< VMGDataGridMap_t, VMGDataGridMap_t> >::const_iterator potentialiter
          = longrangeresults.Result_perIndexSet_LongRange_Grid.find(index);
      if (potentialiter == longrangeresults.Result_perIndexSet_LongRange_Grid.end()) {
        ELOG(1, "appendToHomologyFile() - cannot find index " << (*index)
            << " in FragmentResults.");
        return false;
      }
      // add e+n potentials
      value.potential_distribution =
          boost::fusion::at_key<VMGDataFused::both_sampled_potential>(potentialiter->second.second); // contributions
//      // and re-average to zero (integral is times volume_element which we don't need here)
//      const double sum =
//          std::accumulate(
//              value.potential_distribution.sampled_grid.begin(),
//              value.potential_distribution.sampled_grid.end(),
//              0.);
//      const double offset = sum/(double)value.potential_distribution.sampled_grid.size();
//      for (SamplingGrid::sampledvalues_t::iterator iter = value.potential_distribution.sampled_grid.begin();
//          iter != value.potential_distribution.sampled_grid.end();
//          ++iter)
//        *iter -= offset;
#else
      ELOG(2, "Long-range information in homology desired but long-range analysis capability not compiled in.");
#endif
    }
    values.insert( std::make_pair( graph, value) );
  }
  homology_container.insert(values);

  if (DoLog(2)) {
    LOG(2, "DEBUG: Listing all present atomic ids ...");
    std::stringstream output;
    for (World::AtomIterator iter = World::getInstance().getAtomIter();
        iter != World::getInstance().atomEnd(); ++iter)
      output << (*iter)->getId() << " ";
    LOG(2, "DEBUG: { " << output.str() << "} .");
  }

  // for debugging: print container
  if (DoLog(2)) {
    LOG(2, "DEBUG: Listing all present homologies ...");
    for (HomologyContainer::container_t::const_iterator iter =
        homology_container.begin(); iter != homology_container.end(); ++iter) {
      std::stringstream output;
      output << "DEBUG: graph " << iter->first
          << " has Fragment " << iter->second.fragment
          << ", associated energy " << iter->second.energy;
      if (iter->second.containsGrids)
#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
        output  << ", and sampled grid integral " << iter->second.charge_distribution.integral();
#else
        output << ", and there are sampled grids but capability not compiled in";
#endif
      output << ".";
      LOG(2, output.str());
    }
  }

  return true;
}

// this it taken from
// http://stackoverflow.com/questions/2291802/is-there-a-c-iterator-that-can-iterate-over-a-file-line-by-line
namespace detail
{
  /** Extend the string class by a friend function.
   *
   */
  class Line : public std::string
  {
    friend std::istream & operator>>(std::istream & is, Line & line)
    {
      return std::getline(is, line);
    }
  };
}

/** Parse the given stream line-by-line, passing each to \a dest.
 *
 * \param is stream to parse line-wise
 * \param dest output iterator
 */
template<class OutIt>
void read_lines(std::istream& is, OutIt dest)
{
    typedef std::istream_iterator<detail::Line> InIt;
    std::copy(InIt(is), InIt(), dest);
}


/** Determines the largest cycle in the container and returns its size.
 *
 * \param cycles set of cycles
 * \return size if largest cycle
 */
size_t getMaxCycleSize(const KeySetsContainer &cycles)
{
  // gather cycle sizes
  std::vector<size_t> cyclesizes(cycles.KeySets.size());
  std::transform(
      cycles.KeySets.begin(), cycles.KeySets.end(),
      cyclesizes.begin(),
      boost::bind(&KeySetsContainer::IntVector::size, boost::lambda::_1)
  );
  // get maximum
  std::vector<size_t>::const_iterator maximum_size =
      std::max_element(cyclesizes.begin(), cyclesizes.end());
  if (maximum_size != cyclesizes.end())
    return *maximum_size;
  else
    return 0;
}

void calculateCycleFullContribution(
    const std::map<JobId_t, MPQCData> &shortrangedata,
    const KeySetsContainer &keysets,
    const KeySetsContainer &forcekeysets,
    const KeySetsContainer &cycles,
    const FragmentationShortRangeResults &shortrangeresults)
{
  // copy the shortrangeresults such that private MaxLevel is set in
  // FragmentationShortRangeResults
  FragmentationShortRangeResults cycleresults(shortrangeresults);
  // get largest size
  const size_t maximum_size = getMaxCycleSize(cycles);

  /// The idea here is that (Orthogonal)Summation will place a result
  /// consisting of level 1,2, and 3 fragment and a level 6 ring nonetheless
  /// in level 6. If we want to have this result already at level 3, we
  /// have to specifically inhibit all fragments from later levels but the
  /// cycles and then pick the result from the last level and placing it at
  /// the desired one

  // loop from level 1 to max ring size and gather contributions
  for (size_t level = 1; level <= maximum_size; ++level) {
    // create ValueMask for this level by stepping through each keyset and checking size
    std::vector<bool> localValueMask(shortrangedata.size(), false);
    size_t index=0;
    // TODO: if only KeySetsContainer was usable as a compliant STL container, might be able to use set_difference or alike.
    KeySetsContainer::ArrayOfIntVectors::const_iterator keysetsiter = keysets.KeySets.begin();
    KeySetsContainer::ArrayOfIntVectors::const_iterator cyclesiter = cycles.KeySets.begin();
    for (; (keysetsiter != keysets.KeySets.end()) && (cyclesiter != cycles.KeySets.end());) {
      if (cyclesiter->size() > keysetsiter->size()) {
        // add if not greater than level in size
        if ((*keysetsiter).size() <= level)
          localValueMask[index] = true;
        ++keysetsiter;
        ++index;
      } else if (cyclesiter->size() < keysetsiter->size()) {
        ++cyclesiter;
      } else { // both sets have same size
        if (*cyclesiter > *keysetsiter) {
          // add if not greater than level in size
          if ((*keysetsiter).size() <= level)
            localValueMask[index] = true;
          ++keysetsiter;
          ++index;
        } else if (*cyclesiter < *keysetsiter) {
          ++cyclesiter;
        } else {
          // also always add all cycles
          localValueMask[index] = true;
          ++cyclesiter;
          ++keysetsiter;
          ++index;
        }
      }
    }
    // activate rest if desired by level
    for (; keysetsiter != keysets.KeySets.end(); ++keysetsiter) {
      if ((*keysetsiter).size() <= level)
        localValueMask[index] = true;
      ++index;
    }
    LOG(2, "DEBUG: ValueMask for cycle correction at level " << level << " is "
        << localValueMask << ".");
    // create FragmentationShortRangeResults
    FragmentationShortRangeResults localresults(shortrangedata, keysets, forcekeysets, localValueMask);
    // and perform summation
    localresults(shortrangedata);
    // finally, extract the corrections from last level
    cycleresults.Result_Energy_fused[level-1] =
        localresults.Result_Energy_fused.back();
    cycleresults.Result_Time_fused[level-1] =
        localresults.Result_Time_fused.back();
    cycleresults.Result_Force_fused[level-1] =
        localresults.Result_Force_fused.back();
  }
  printReceivedCycleResults(cycleresults);
}

static void AddForces(
    const IndexedVectors::indexedvectors_t &_forces,
    const bool _IsAngstroem)
{
  for(IndexedVectors::indexedvectors_t::const_iterator iter = _forces.begin();
      iter != _forces.end(); ++iter) {
    const IndexedVectors::index_t &index = iter->first;
    const IndexedVectors::vector_t &forcevector = iter->second;
    ASSERT( forcevector.size() == NDIM,
        "printReceivedShortResults() - obtained force vector has incorrect dimension.");
    // note that mpqc calculates a gradient, hence force pointing into opposite direction
    // we have to mind different units here: MPQC has a_o, while we may have angstroem
    Vector ForceVector(-forcevector[0], -forcevector[1], -forcevector[2]);
    if (_IsAngstroem)
      for (size_t i=0;i<NDIM;++i)
        ForceVector[i] *= AtomicLengthToAngstroem;
    atom *_atom = World::getInstance().getAtom(AtomById(index));
    if(_atom != NULL)
      _atom->setAtomicForce(_atom->getAtomicForce() + ForceVector);
    else
      ELOG(2, "Could not find atom to id " << index << ".");
  }
}

ActionState::ptr FragmentationAnalyseFragmentationResultsAction::performCall() {
  bool status=true;

  /// get data and keysets from ResultsContainer
  FragmentationResultContainer& container = FragmentationResultContainer::getInstance();
  const std::map<JobId_t, MPQCData> &shortrangedata = container.getShortRangeResults();
  const KeySetsContainer &keysets = container.getKeySets();
  const KeySetsContainer &forcekeysets = container.getForceKeySets();
  const bool DoLongrange = container.areFullRangeResultsPresent();
  const bool IsAngstroem = true;

  if (keysets.KeySets.empty()) {
    STATUS("There are no results in the container.");
    return Action::failure;
  }

  /// calculate normal contributions with (if present) cycles coming at their
  /// respective bond order.
  std::vector<bool> ValueMask(shortrangedata.size(), true);
  FragmentationShortRangeResults shortrangeresults(shortrangedata, keysets, forcekeysets, ValueMask);
  shortrangeresults(shortrangedata);
  printReceivedShortResults(shortrangeresults);
  printReceivedShortResultsPerIndex(shortrangeresults);
  // add summed results to container
  container.addShortRangeSummedResults(shortrangeresults.getSummedShortRangeResults());

  /// now do we need to calculate the cycle contribution
  // check whether there are cycles in container or else in file
  KeySetsContainer cycles = container.getCycles();
  if (cycles.KeySets.empty()) {
    // parse from file if cycles is empty
    boost::filesystem::path filename(
        params.prefix.get() + std::string(CYCLEKEYSETFILE));
    if (boost::filesystem::exists(filename)) {
      LOG(1, "INFO: Parsing cycles file " << filename.string() << ".");
      // parse file line by line
      std::ifstream File;
      File.open(filename.string().c_str());
      typedef std::istream_iterator<detail::Line> InIt;
      for (InIt iter = InIt(File); iter != InIt(); ++iter) {
        KeySetsContainer::IntVector cycle;
        std::stringstream line(*iter);
        while (line.good()) {
          int id;
          line >> id >> ws;
          cycle.push_back(id);
        }
        if (!cycle.empty()) {
          LOG(2, "DEBUG: Adding cycle " << cycle << ".");
          cycles.insert( cycle, cycle.size());
        }
      }
      File.close();
    } else {
      LOG(1, "INFO: Cycles file not found at " << filename.string() << ".");
    }
  }

  // calculate energy if cycles are calculated fully at each level already
  if (!cycles.KeySets.empty()) {
    calculateCycleFullContribution(
        shortrangedata,
        keysets,
        forcekeysets,
        cycles,
        shortrangeresults);
  }

  // adding obtained forces
  if ( const_cast<const World &>(World::getInstance()).getAllAtoms().size() != 0) {
    const IndexedVectors::indexedvectors_t shortrange_forces =
        boost::fusion::at_key<MPQCDataFused::forces>(
            shortrangeresults.Result_Force_fused.back()
            ).getVectors();
    AddForces(shortrange_forces,IsAngstroem);
  } else {
    LOG(1, "INFO: Full molecule not loaded, hence will not add forces to atoms.");
  }

#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
  if (DoLongrange) {
    if ( const_cast<const World &>(World::getInstance()).getAllAtoms().size() == 0) {
      STATUS("Please load the full molecule into std::map<JobId_t, VMGData> longrangeData the world before starting this action.");
      return Action::failure;
    }

    std::map<JobId_t, VMGData> longrangeData = container.getLongRangeResults();

    ASSERT( !longrangeData.empty(),
        "FragmentationAnalyseFragmentationResultsAction::performCall() - longrangeData is empty?");
    // get required level from longrange full potential solution
    const int level = longrangeData.rbegin()->second.sampled_potential.level;
    SamplingGridProperties domain(ExportGraph_ToJobs::getDomainGrid(level));
    SamplingGrid zero_globalgrid(domain);
    FragmentationChargeDensity summedChargeDensity(shortrangedata);
    summedChargeDensity(shortrangedata, keysets, zero_globalgrid);
    const std::vector<SamplingGrid> full_sample = summedChargeDensity.getFullSampledGrid();

    // remove full solution corresponding to full_sample from map (must be highest ids), has to be treated extra
    std::map<JobId_t, VMGData>::iterator iter = longrangeData.end();
    std::advance(iter, -full_sample.size());
    std::map<JobId_t, VMGData>::iterator remove_iter = iter;
    std::vector<VMGData> fullsolutionData;
    for (; iter != longrangeData.end(); ++iter)
      fullsolutionData.push_back(iter->second);
    if (longrangeData.size() > 1) // when there's just a single fragment, it corresponds to full solution
      longrangeData.erase(remove_iter, longrangeData.end());

    // Final phase: sum up and print result
    IndexedVectors::indices_t implicit_indices;
    if (params.UseImplicitCharges.get()) {
      // place all in implicit charges that are not selected but contained in ParticleRegistry
      const World &world = const_cast<const World &>(World::getInstance());
      const ParticleRegistry &registry = const_cast<const ParticleRegistry &>(ParticleRegistry::getInstance());
      const World::ConstAtomComposite &atoms = world.getAllAtoms();
      for (World::ConstAtomComposite::const_iterator iter = atoms.begin();
          iter != atoms.end(); ++iter) {
        const atomId_t atomid = (*iter)->getId();
        if (!world.isAtomSelected(atomid)) {
          const std::string &symbol = (*iter)->getElement().getSymbol();
          if (registry.isPresentByName(symbol))
            implicit_indices.push_back(atomid);
        }
      }
      LOG(2, "INFO: We added " << implicit_indices.size() << " indices due to implicit charges.");
    }

    FragmentationLongRangeResults longrangeresults(
        shortrangedata, longrangeData, keysets, forcekeysets);
    {
      SamplingGrid zero_globalgrid(domain);
      longrangeresults(
          shortrangedata,
          longrangeData,
          fullsolutionData,
          full_sample,
          zero_globalgrid,
          implicit_indices);
      printReceivedFullResults(longrangeresults);
    }

    // add long-range forces
    if ( const_cast<const World &>(World::getInstance()).getAllAtoms().size() != 0) {
      const IndexedVectors::indexedvectors_t longrange_forces =
          boost::fusion::at_key<VMGDataFused::forces_longrange>(
              longrangeresults.Result_ForcesLongRangeIntegrated_fused.back()
              ).getVectors();
      AddForces(longrange_forces,IsAngstroem);
    } else {
      LOG(1, "INFO: Full molecule not loaded, hence will not add forces to atoms.");
    }

    // append all keysets to homology file
    status = appendToHomologies(shortrangeresults, longrangeresults, params.DoStoreGrids.get());
  } else {
    // append all keysets to homology file with short-range info only (without grids)
    std::map<JobId_t, VMGData> longrangeData;
    FragmentationLongRangeResults longrangeresults(
        shortrangedata, longrangeData, keysets, forcekeysets);
    status = appendToHomologies(shortrangeresults, longrangeresults, false);
  }
#else
  if (DoLongrange) {
    ELOG(2, "File contains long-range information but long-range analysis capability not compiled in.");
  }

  // append all keysets to homology file with short-range info only (without grids)
  status = appendToHomologies(shortrangeresults, false);
#endif

  // we no longer clear the container
//  container.clear();

  if (status)
    return Action::success;
  else {
    STATUS("AnalyseFragmentResultsAction failed: invalid results, failed to append to homologies.");
    return Action::failure;
  }
}

ActionState::ptr FragmentationAnalyseFragmentationResultsAction::performUndo(ActionState::ptr _state) {
  return Action::success;
}

ActionState::ptr FragmentationAnalyseFragmentationResultsAction::performRedo(ActionState::ptr _state){
  return Action::success;
}

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

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