Changeset 3690e4 for src/Fragmentation

Timestamp:

Feb 2, 2016, 5:50:29 PM (9 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

9ae11c

Parents:

d1831e (diff), 62d092 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'Enhancing_Interdistance' into Candidate_v1.5.1

Location:

src/Fragmentation

Files:

• Automation/ResultContainer_impl.hpp (modified) (1 diff)
• Automation/VMGFragmentController.cpp (modified) (3 diffs)
• Automation/VMGFragmentController.hpp (modified) (2 diffs)
• Exporters/ExportGraph_ToJobs.cpp (modified) (1 diff)
• Interfragmenter.cpp (modified) (7 diffs)
• Interfragmenter.hpp (modified) (3 diffs)
• Makefile.am (modified) (1 diff)
• Summation/Containers/FragmentationLongRangeResults.cpp (modified) (5 diffs)
• Summation/Containers/FragmentationLongRangeResults.hpp (modified) (2 diffs)
• Summation/Containers/FragmentationShortRangeResults.cpp (modified) (1 diff)
• Summation/Containers/MPQCData.hpp (modified) (4 diffs)
• Summation/Containers/MPQCDataFused.hpp (modified) (2 diffs)
• Summation/Containers/VMGData.cpp (modified) (2 diffs)
• Summation/Containers/VMGData.hpp (modified) (3 diffs)
• Summation/Containers/VMGDataFused.hpp (modified) (3 diffs)
• Summation/Containers/VMGDataMap.hpp (modified) (4 diffs)
• Summation/Containers/VMGData_printKeyNames.hpp (modified) (2 diffs)
• Summation/Converter/DataConverter.hpp (modified) (3 diffs)
• Summation/Makefile.am (modified) (1 diff)
• Summation/SetValues/FragmentForces.cpp (added)
• Summation/SetValues/FragmentForces.hpp (added)
• Summation/SetValues/Makefile.am (modified) (2 diffs)

src/Fragmentation/Automation/ResultContainer_impl.hpp

@@ -38 +38 @@
     boost::archive::text_iarchive ia(inputstream);
     ia >> extractedData;
-    LOG(1, "INFO: extracted data is " << extractedData << ".");
+    LOG(1, "INFO: extracted data of #" << (*iter)->getId() << " is " << extractedData << ".");
     fragmentData.push_back(extractedData);
   }

src/Fragmentation/Automation/VMGFragmentController.cpp

@@ -76 +76 @@
     const MPQCData::DoValenceOnly_t _DoValenceOnly,
     const bool _DoPrintDebug,
-    const bool _OpenBoundaryConditions)
+    const bool _OpenBoundaryConditions,
+    const bool _DoSmearCharges)
 {
   std::vector<FragmentJob::ptr> jobs;
@@ -99 +100 @@
             _SampleParticles == DoSampleParticles,
             _DoPrintDebug,
-            _OpenBoundaryConditions) );
+            _OpenBoundaryConditions,
+            _DoSmearCharges) );
     jobs.push_back(testJob);
   }
@@ -147 +149 @@
             _SampleParticles == DoSampleParticles,
             _DoPrintDebug,
-            _OpenBoundaryConditions) );
+            _OpenBoundaryConditions,
+            _DoSmearCharges) );
     jobs.push_back(testJob);
   }

src/Fragmentation/Automation/VMGFragmentController.hpp

@@ -63 +63 @@
    * \param _OpenBoundaryConditions whether we have open (true) or periodic (false)
    *        boundary conditions
+   * \param _DoSmearCharges whether to smear out electronic charge distributions with bsplines or not
    */
   bool createLongRangeJobs(
@@ -73 +74 @@
       const MPQCData::DoValenceOnly_t _DoValenceOnly,
       const bool _DoPrintDebug,
-      const bool _OpenBoundaryConditions = false);
+      const bool _OpenBoundaryConditions = false,
+      const bool _DoSmearCharges = false);
 
   void waitforResults(const size_t NoExpectedResults)

src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

@@ -80 +80 @@
   RealSpaceMatrix M = World::getInstance().getDomain().getM();
   M *= 1./AtomicLengthToAngstroem;  // scale to atomic length units
-  const double size = M.at(0,0);
+  const double size = std::max( std::max(M.at(0,0), M.at(1,1)), M.at(2,2));
   double end[NDIM] = { size, size, size };
   const ParserTypes jobtype =

src/Fragmentation/Interfragmenter.cpp

@@ -37 +37 @@
 #include "Interfragmenter.hpp"
 
-#include <list>
-#include <map>
-
 #include "CodePatterns/Assert.hpp"
 #include "CodePatterns/Log.hpp"
@@ -45 +42 @@
 #include "LinearAlgebra/Vector.hpp"
 
-#include "AtomIdSet.hpp"
+#include "Descriptors/AtomDescriptor.hpp"
 #include "Element/element.hpp"
 #include "Fragmentation/Graph.hpp"
@@ -53 +50 @@
 #include "World.hpp"
 
-void Interfragmenter::operator()(
-    size_t MaxOrder,
-    double Rcut,
-    const enum HydrogenTreatment treatment)
+
+Interfragmenter::atomkeyset_t Interfragmenter::getAtomKeySetMap(
+    size_t _MaxOrder
+    ) const
 {
   /// create a map of atom to keyset (below equal MaxOrder)
-  typedef std::list<const KeySet *> keysets_t;
-  typedef std::map<const atom *, keysets_t > atomkeyset_t;
   atomkeyset_t atomkeyset;
   LOG(1, "INFO: Placing all atoms and their keysets into a map.");
@@ -69 +64 @@
     const AtomIdSet atoms(keyset);
     const size_t atoms_size = atoms.getAtomIds().size();
-    if ((atoms_size > MaxOrder) || (atoms_size == 0))
+    if ((atoms_size > _MaxOrder) || (atoms_size == 0))
       continue;
     for (AtomIdSet::const_iterator atomiter = atoms.begin();
@@ -87 +82 @@
   LOG(2, "DEBUG: There are " << atomkeyset.size() << " entries in lookup.");
 
+  return atomkeyset;
+}
+
+static Vector getAtomIdSetCenter(
+    const AtomIdSet &_atoms)
+{
+  const molecule * const _mol = (*_atoms.begin())->getMolecule();
+  const size_t atoms_size = _atoms.getAtomIds().size();
+  Vector center;
+  for (AtomIdSet::const_iterator iter = _atoms.begin();
+      iter != _atoms.end(); ++iter) {
+    center += (*iter)->getPosition();
+    ASSERT ( _mol == (*iter)->getMolecule(),
+        "getAtomIdSetCenter() - ids in same keyset belong to different molecule.");
+  }
+  center *= 1./(double)atoms_size;
+
+  return center;
+}
+
+Interfragmenter::candidates_t Interfragmenter::getNeighborsOutsideMolecule(
+    const AtomIdSet &_atoms,
+    double _Rcut,
+    const enum HydrogenTreatment _treatment) const
+{
+  /// go through linked cell and get all neighboring atoms up to Rcut
+  const LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(_Rcut);
+  const Vector center = getAtomIdSetCenter(_atoms);
+  const LinkedCell::LinkedList neighbors = view.getAllNeighbors(_Rcut, center);
+  LOG(4, "DEBUG: Obtained " << neighbors.size() << " neighbors in distance of "
+      << _Rcut << " around " << center << ".");
+
+  /// remove all atoms that belong to same molecule as does one of the
+  /// fragment's atoms
+  const molecule * const _mol = (*_atoms.begin())->getMolecule();
+  candidates_t candidates;
+  candidates.reserve(neighbors.size());
+  for (LinkedCell::LinkedList::const_iterator iter = neighbors.begin();
+      iter != neighbors.end(); ++iter) {
+    const atom * const _atom = static_cast<const atom * const >(*iter);
+    ASSERT( _atom != NULL,
+        "Interfragmenter::getNeighborsOutsideMolecule() - a neighbor is not actually an atom?");
+    if ((_atom->getMolecule() != _mol)
+        && (_atom->getPosition().DistanceSquared(center) < _Rcut*_Rcut)
+        && ((_treatment == IncludeHydrogen) || (_atom->getType()->getAtomicNumber() != 1))) {
+      candidates.push_back(_atom);
+    }
+  }
+  LOG(3, "DEBUG: There remain " << candidates.size() << " candidates.");
+
+  return candidates;
+}
+
+Interfragmenter::atomkeyset_t Interfragmenter::getCandidatesSpecificKeySetMap(
+    const candidates_t &_candidates,
+    const atomkeyset_t &_atomkeyset) const
+{
+  atomkeyset_t fragmentmap;
+  for (candidates_t::const_iterator candidateiter = _candidates.begin();
+      candidateiter != _candidates.end(); ++candidateiter) {
+    const atom * _atom = *candidateiter;
+    atomkeyset_t::const_iterator iter = _atomkeyset.find(_atom);
+    ASSERT( iter != _atomkeyset.end(),
+        "Interfragmenter::getAtomSpecificKeySetMap() - could not find atom "
+        +toString(_atom->getId())+" in lookup.");
+    fragmentmap.insert( std::make_pair( _atom, iter->second) );
+  }
+  LOG(4, "DEBUG: Copied part of lookup map contains " << fragmentmap.size() << " keys.");
+
+  return fragmentmap;
+}
+
+void Interfragmenter::combineFragments(
+    const size_t _MaxOrder,
+    const candidates_t &_candidates,
+    atomkeyset_t &_fragmentmap,
+    const KeySet &_keyset,
+    Graph &_InterFragments,
+    int &_counter)
+{
+  for (candidates_t::const_iterator candidateiter = _candidates.begin();
+      candidateiter != _candidates.end(); ++candidateiter) {
+    const atom *_atom = *candidateiter;
+    LOG(3, "DEBUG: Current candidate is " << *_atom << ".");
+    atomkeyset_t::iterator finditer = _fragmentmap.find(_atom);
+    ASSERT( finditer != _fragmentmap.end(),
+        "Interfragmenter::combineFragments() - could not find atom "
+        +toString(_atom->getId())+" in fragment specific lookup.");
+    keysets_t &othersets = finditer->second;
+    ASSERT( !othersets.empty(),
+        "Interfragmenter::combineFragments() - keysets to "+toString(_atom->getId())+
+        "is empty.");
+    keysets_t::iterator otheriter = othersets.begin();
+    while (otheriter != othersets.end()) {
+      const KeySet &otherset = **otheriter;
+      LOG(3, "DEBUG: Current keyset is " << otherset << ".");
+      // only add them one way round and not the other
+      if (otherset < _keyset) {
+        ++otheriter;
+        continue;
+      }
+      // only add if combined they don't exceed the desired maxorder
+      if (otherset.size() + _keyset.size() > _MaxOrder) {
+        LOG(3, "INFO: Rejecting " << otherset << " as in sum their orders exceed " << _MaxOrder);
+        ++otheriter;
+        continue;
+      }
+      KeySet newset(otherset);
+      newset.insert(_keyset.begin(), _keyset.end());
+      LOG(3, "DEBUG: Inserting new combined set " << newset << ".");
+      _InterFragments.insert( std::make_pair(newset, std::make_pair(++_counter, 1.)));
+      // finally, remove the set such that no other combination exists
+      otheriter = othersets.erase(otheriter);
+    }
+  }
+}
+
+void Interfragmenter::operator()(
+    const size_t MaxOrder,
+    const double Rcut,
+    const enum HydrogenTreatment treatment)
+{
+  atomkeyset_t atomkeyset = getAtomKeySetMap(MaxOrder);
+
   Graph InterFragments;
   int counter = TotalGraph.size();
@@ -101 +220 @@
       continue;
 
-    /// go through linked cell and get all neighboring atoms up to Rcut
-    Vector center;
-    const molecule *_mol = (*atoms.begin())->getMolecule();
-    for (AtomIdSet::const_iterator iter = atoms.begin();
-        iter != atoms.end(); ++iter) {
-      center += (*iter)->getPosition();
-      ASSERT ( _mol == (*iter)->getMolecule(),
-          "Interfragmenter::operator() - ids in same keyset belong to different molecule.");
-    }
-    center *= 1./(double)atoms_size;
-    LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(Rcut);
-    LinkedCell::LinkedList neighbors = view.getAllNeighbors(Rcut, center);
-    LOG(4, "DEBUG: Obtained " << neighbors.size() << " neighbors in distance of "
-        << Rcut << " around " << center << ".");
-
-    /// remove all atoms that belong to same molecule as does one of the
-    /// fragment's atoms
-    typedef std::vector<const atom *> candidates_t;
-    candidates_t candidates;
-    candidates.reserve(neighbors.size());
-    for (LinkedCell::LinkedList::const_iterator iter = neighbors.begin();
-        iter != neighbors.end(); ++iter) {
-      const atom * const _atom = static_cast<const atom * const >(*iter);
-      ASSERT( _atom != NULL,
-          "Interfragmenter::operator() - a neighbor is not actually an atom?");
-      if ((_atom->getMolecule() != _mol)
-          && (_atom->getPosition().DistanceSquared(center) < Rcut*Rcut)
-          && ((treatment == IncludeHydrogen) || (_atom->getType()->getAtomicNumber() != 1))) {
-        candidates.push_back(_atom);
-      }
-    }
-    LOG(3, "DEBUG: There remain " << candidates.size() << " candidates.");
+    // get neighboring atoms outside the current molecule
+    candidates_t candidates = getNeighborsOutsideMolecule(atoms, Rcut, treatment);
 
     // create a lookup specific to this fragment
-    atomkeyset_t fragmentmap;
-    for (candidates_t::const_iterator candidateiter = candidates.begin();
-        candidateiter != candidates.end(); ++candidateiter) {
-      const atom * _atom = *candidateiter;
-      atomkeyset_t::const_iterator iter = atomkeyset.find(_atom);
-      ASSERT( iter != atomkeyset.end(),
-          "Interfragmenter::operator() - could not find atom "
-          +toString(_atom->getId())+" in lookup.");
-      fragmentmap.insert( std::make_pair( _atom, iter->second) );
-    }
-    LOG(4, "DEBUG: Copied part of lookup map contains " << fragmentmap.size() << " keys.");
+    atomkeyset_t fragmentmap = getCandidatesSpecificKeySetMap(candidates, atomkeyset);
 
     /// combine each remaining fragment with current fragment to a new fragment
-    /// if keyset is less
-    for (candidates_t::const_iterator candidateiter = candidates.begin();
-        candidateiter != candidates.end(); ++candidateiter) {
-      const atom *_atom = *candidateiter;
-      LOG(3, "DEBUG: Current candidate is " << *_atom << ".");
-      atomkeyset_t::iterator finditer = fragmentmap.find(_atom);
-      ASSERT( finditer != fragmentmap.end(),
-          "Interfragmenter::operator() - could not find atom "
-          +toString(_atom->getId())+" in fragment specific lookup.");
-      keysets_t &othersets = finditer->second;
-      keysets_t::iterator otheriter = othersets.begin();
-      while (otheriter != othersets.end()) {
-        const KeySet &otherset = **otheriter;
-        LOG(3, "DEBUG: Current keyset is " << otherset << ".");
-        // only add them one way round and not the other
-        if (otherset < keyset) {
-          ++otheriter;
-          continue;
-        }
-        KeySet newset(otherset);
-        newset.insert(keyset.begin(), keyset.end());
-        LOG(3, "DEBUG: Inserting new combined set " << newset << ".");
-        InterFragments.insert( std::make_pair(newset, std::make_pair(++counter, 1.)));
-        // finally, remove the set such that no other combination exists
-        otheriter = othersets.erase(otheriter);
-      }
-    }
+    /// if keyset is less (to prevent addding same inter-fragment twice)
+    combineFragments(MaxOrder, candidates, fragmentmap, keyset, InterFragments, counter);
   }
 
@@ -181 +235 @@
   TotalGraph.InsertGraph(InterFragments, counter);
 }
+
+double Interfragmenter::findLargestCutoff(
+    const size_t _MaxOrder,
+    const double _upperbound,
+    const enum HydrogenTreatment _treatment) const
+{
+  double Rcut = _upperbound*_upperbound;
+  std::pair<atomId_t, atomId_t> ClosestPair;
+  // place all atoms into LC grid with some upper bound
+  const LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(_upperbound);
+
+  atomkeyset_t atomkeyset = getAtomKeySetMap(_MaxOrder);
+
+  // go through each atom and find closest atom not in the same keyset
+  for (Graph::const_iterator keysetiter = TotalGraph.begin();
+      keysetiter != TotalGraph.end(); ++keysetiter) {
+    const KeySet &keyset = keysetiter->first;
+    LOG(2, "DEBUG: Current keyset is " << keyset);
+    const AtomIdSet atoms(keyset);
+
+    // get neighboring atoms outside the current molecule
+    const candidates_t candidates = getNeighborsOutsideMolecule(atoms, _upperbound, _treatment);
+    const Vector center = getAtomIdSetCenter(atoms);
+
+    for (candidates_t::const_iterator candidateiter = candidates.begin();
+        candidateiter != candidates.end(); ++candidateiter) {
+      atom const * const Walker = *candidateiter;
+      // go through each atom in set and pick minimal distance
+      for (AtomIdSet::const_iterator setiter = atoms.begin(); setiter != atoms.end(); ++setiter) {
+        const double distanceSquared = Walker->getPosition().DistanceSquared(center);
+        // pick the smallest compared to current Rcut if smaller
+        if (distanceSquared < Rcut) {
+          Rcut = distanceSquared;
+          ClosestPair.first = (*setiter)->getId();
+          ClosestPair.second = Walker->getId();
+          LOG(2, "DEBUG: Found new pair " << ClosestPair << " with distance " << sqrt(Rcut));
+        }
+      }
+    }
+  }
+  const double largest_distance = sqrt(Rcut);
+  LOG(1, "INFO: Largest inter-fragment distance to cause no additional fragments: "
+      << largest_distance);
+
+  return largest_distance;
+}

src/Fragmentation/Interfragmenter.hpp

@@ -14 +14 @@
 #endif
 
+#include <list>
+#include <map>
 #include <vector>
 
+#include "AtomIdSet.hpp"
 #include "Fragmentation/HydrogenSaturation_enum.hpp"
 
@@ -50 +53 @@
    */
   void operator()(
-      size_t MaxOrder,
-      double Rcut,
+      const size_t MaxOrder,
+      const double Rcut,
       const enum HydrogenTreatment treatment);
+
+  /** This finds the largest cut off distance (\a Rcut) such that when running
+   * operator() no additional inter-fragments would be produced.
+   *
+   * \param MaxOrder maximum order for fragments to interrelate
+   * \param _upperbound upper bound on \a Rcut above which we do not look
+   * \param treatment whether hydrogens are treated specially or not
+   * \return largest cutoff distance to cause no additional inter-fragments
+   */
+  double findLargestCutoff(
+      const size_t _MaxOrder,
+      const double _upperbound,
+      const enum HydrogenTreatment _treatment) const;
 
 private:
@@ -62 +78 @@
   std::vector<atom *> getAtomsFromKeySet(const KeySet &keyset) const;
 
+  typedef std::list<const KeySet *> keysets_t;
+  typedef std::map<const atom *, keysets_t > atomkeyset_t;
+
+  /** Helper function to create a map of atoms and their fragments/keysets.
+   *
+   * \param _MaxOrder maximum order
+   * \return map with atoms as a keys and fragments as values
+   */
+  atomkeyset_t getAtomKeySetMap(size_t _MaxOrder) const;
+
+  typedef std::vector<const atom *> candidates_t;
+
+  /** Helper function to get all atoms around a specific keyset not contained in
+   * the same molecule.
+   *
+   * \param _atoms all atoms of a fragment
+   * \param _Rcut desired distance cutoff
+   * \param _treatment whether hydrogens are treated special or not
+   */
+  candidates_t getNeighborsOutsideMolecule(
+      const AtomIdSet &_atoms,
+      const double _Rcut,
+      const enum HydrogenTreatment _treatment) const;
+
+  /** Helper function to return a fragment/KeySet map specific to all candidates.
+   *
+   * \param _candidates all neighboring atoms around keyset
+   * \param _atomkeyset map with all atoms and the KeySets they are contained in
+   * \return specific fragment map
+   */
+  atomkeyset_t getCandidatesSpecificKeySetMap(
+      const candidates_t &_candidates,
+      const atomkeyset_t &_atomkeyset) const;
+
+  /** For a given set of candidates atoms in \a _candidates and a \a _keyset
+   * we combine each fragment from either atom and place it into internal
+   * Graph.
+   *
+   * \param _MaxOrder maximum order
+   * \param _candidates all atoms neighboring the current out outside of its molecule
+   * \param _fragmentmap all keysets related to this atom
+   * \param _keyset current keyset (used as base for creating inter-fragments)
+   * \param _InterFragments container for all created inter-fragments
+   * \param _counter counts added fragments
+   */
+  void combineFragments(
+      const size_t _MaxOrder,
+      const candidates_t &_candidates,
+      atomkeyset_t &_fragmentmap,
+      const KeySet &_keyset,
+      Graph &_InterFragments,
+      int &_counter);
 
 private:

src/Fragmentation/Makefile.am

@@ -128 +128 @@
 libMolecuilderFragmentation_KeysetsContainer_libincludedir = $(libdir)/MoleCuilder/include
 nodist_libMolecuilderFragmentation_KeysetsContainer_libinclude_HEADERS = $(top_builddir)/libmolecuilder_config.h
-
-## Install the generated pkg-config file (.pc) into the expected location for
-## architecture-dependent package configuration information.  Occasionally,
-## pkg-config files are also used for architecture-independent data packages,
-## in which case the correct install location would be $(datadir)/pkgconfig.
-#pkgconfigdir = $(libdir)/pkgconfig
-#pkgconfig_DATA = $(top_builddir)/MoleCuilder.pc

src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp

@@ -65 +65 @@
     const KeySetsContainer& _ForceKeySet) :
     KeySet(_KeySet),
-    ForceKeySet(_ForceKeySet)
+    ForceKeySet(_ForceKeySet),
+    hasForces((!longrangeData.empty()) && (longrangeData.begin()->second.hasForces))
 {
   initLookups(fragmentData, longrangeData);
@@ -126 +127 @@
         Result_LongRange_fused, Result_perIndexSet_LongRange);
 
+    IndexedVectors::indices_t fullindices;
+    if (hasLongRangeForces()) {
+      // force has extra data converter (this is similar to MPQCData's forces
+      std::map<JobId_t, VMGDataForceMap_t> VMGData_Force_fused;
+      convertDatatoForceMap<VMGData, VMGDataForceMap_t, VMGDataFused>(
+          longrangeData, ForceKeySet, VMGData_Force_fused);
+      Result_ForceLongRange_fused.resize(MaxLevel); // we need the results of correct size already
+      AllLevelOrthogonalSummator<VMGDataForceMap_t> forceSummer(
+                  subsetmap,
+                  VMGData_Force_fused,
+                  container->getContainer(),
+                  VMGMatrixNrLookup,
+                  Result_ForceLongRange_fused,
+                  Result_perIndexSet_LongRange_Force);
+      boost::mpl::for_each<VMGDataForceVector_t>(boost::ref(forceSummer));
+      // build full force index set
+      KeySetsContainer::ArrayOfIntVectors::const_iterator arrayiter = ForceKeySet.KeySets.begin();
+      std::set<IndexedVectors::index_t> sorted_indices;
+      for (;arrayiter != ForceKeySet.KeySets.end(); ++arrayiter) {
+        sorted_indices.insert(arrayiter->begin(), arrayiter->end());
+      }
+      sorted_indices.erase(-1);
+      fullindices.insert(fullindices.begin(), sorted_indices.begin(), sorted_indices.end());
+    }
+
+    // then sum up
+    OrthogonalSumUpPerLevel<VMGDataGridMap_t, VMGData, VMGDataGridVector_t>(
+        longrangeData, VMGMatrixNrLookup, container, subsetmap,
+        Result_GridLongRange_fused, Result_perIndexSet_LongRange_Grid);
+
     Result_LongRangeIntegrated_fused.reserve(MaxLevel);
     // NOTE: potential for level 1 is wrongly calculated within a molecule
@@ -146 +177 @@
       // the full solution.
       const SamplingGrid &short_range_correction =
-          boost::fusion::at_key<VMGDataFused::sampled_potential>(Result_LongRange_fused[level-1]);
+          boost::fusion::at_key<VMGDataFused::sampled_potential>(Result_GridLongRange_fused[level-1]);
       double electron_short_range_energy = short_range_correction.integral();
       full_sample_solution -= short_range_correction;
@@ -202 +233 @@
           + boost::fusion::at_key<VMGDataFused::nuclei_shortrange>(instance);
       Result_LongRangeIntegrated_fused.push_back(instance);
+
+      if (hasLongRangeForces()) {
+        VMGDataLongRangeForceMap_t forceinstance;
+        IndexedVectors fullforces( fullindices, fullsolutionData[level-1].forces);
+        IndexedVectors longrangeforces =
+            boost::fusion::at_key<VMGDataFused::forces>(Result_ForceLongRange_fused[level-1]);
+        boost::fusion::at_key<VMGDataFused::forces_shortrange>(forceinstance) =
+            fullforces;
+        fullforces -= longrangeforces;
+        boost::fusion::at_key<VMGDataFused::forces_longrange>(forceinstance) =
+            fullforces;
+        Result_ForcesLongRangeIntegrated_fused.push_back(forceinstance);
+      }
     }
 //    {
@@ -207 +251 @@
 //      SamplingGrid full_sample_solution = fullsolutionData.back().sampled_potential;
 //      const SamplingGrid &short_range_correction =
-//          boost::fusion::at_key<VMGDataFused::sampled_potential>(Result_LongRange_fused.back()).getFullContribution();
+//          boost::fusion::at_key<VMGDataFused::sampled_potential>(Result_GridLongRange_fused.back()).getFullContribution();
 //      full_sample_solution -= short_range_correction;
 //      // multiply element-wise with charge distribution

src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp

@@ -87 +87 @@
   }
 
+  const bool hasLongRangeForces() const
+  { return hasForces; }
+
 private:
   void initLookups(
@@ -100 +103 @@
   IndexSetContainer::ptr container;
   SubsetMap::ptr subsetmap;
+  //!> indicates whether the longrange results contain VMG forces data or not
+  bool hasForces;
 
 public:
   //!> results per level of summed up sampled grid charge
   std::vector<MPQCDataGridMap_t> Result_Grid_fused;
-  //!> results per level of summed up long range potential grids and energy
+  //!> results per level of summed up long range energies
   std::vector<VMGDataMap_t> Result_LongRange_fused;
+  //!> results per level of summed up long range forces
+  std::vector<VMGDataForceMap_t> Result_ForceLongRange_fused;
+  //!> results per level of summed up long range potential grids
+  std::vector<VMGDataGridMap_t> Result_GridLongRange_fused;
   //!> results per level of summed up long range true energy
   std::vector<VMGDataLongRangeMap_t> Result_LongRangeIntegrated_fused;
+  //!> results per level of summed up long range true forces
+  std::vector<VMGDataLongRangeForceMap_t> Result_ForcesLongRangeIntegrated_fused;
 
   //!> results per IndexSet of summed up sampled grid charge
   std::map<IndexSet::ptr, std::pair<MPQCDataGridMap_t,MPQCDataGridMap_t> > Result_perIndexSet_Grid;
-  //!> results per IndexSet of summed up long range potential grids and energy
+  //!> results per IndexSet of summed up long range energies
   std::map<IndexSet::ptr, std::pair<VMGDataMap_t, VMGDataMap_t> > Result_perIndexSet_LongRange;
+  //!> results per IndexSet of summed up long range forces
+  std::map<IndexSet::ptr, std::pair<VMGDataForceMap_t, VMGDataForceMap_t> > Result_perIndexSet_LongRange_Force;
+  //!> results per IndexSet of summed up long range potential grids
+  std::map<IndexSet::ptr, std::pair<VMGDataGridMap_t, VMGDataGridMap_t> > Result_perIndexSet_LongRange_Grid;
   // we don't need the map pendant for Result_LongRangeIntegrated_fused, as this
   // quantity makes sense only level-wise

src/Fragmentation/Summation/Containers/FragmentationShortRangeResults.cpp

@@ -109 +109 @@
     // force has extra data converter
     std::map<JobId_t, MPQCDataForceMap_t> MPQCData_Force_fused;
-    convertMPQCDatatoForceMap(fragmentData, ForceKeySet, MPQCData_Force_fused);
+    convertDatatoForceMap<MPQCData, MPQCDataForceMap_t, MPQCDataFused>(fragmentData, ForceKeySet, MPQCData_Force_fused);
     Result_Force_fused.resize(MaxLevel); // we need the results of correct size already
     AllLevelOrthogonalSummator<MPQCDataForceMap_t> forceSummer(

src/Fragmentation/Summation/Containers/MPQCData.hpp

@@ -22 +22 @@
 
 #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
+#include "Fragmentation/Summation/SetValues/FragmentForces.hpp"
 
 class MPQCCommandJob;
@@ -77 +78 @@
 
   /// Forces
-  typedef std::vector<double> vector_type;
-  std::vector< vector_type  > forces;
+  FragmentForces forces;
 
   //!> whether to actually sample the density
@@ -131 +131 @@
     ar & energies.hcore;
     ar & energies.eigenvalues;
-    ar & forces;
+    if (version < 2) {
+      FragmentForces::forces_t justforces(forces);
+      ar & justforces;
+      dynamic_cast<FragmentForces::forces_t &>(forces) = justforces;
+    } else
+      ar & forces;
     ar & sampled_grid;
     ar & DoLongrange;
@@ -147 +152 @@
 };
 
-BOOST_CLASS_VERSION(MPQCData, 1)
+BOOST_CLASS_VERSION(MPQCData, 2)
 
 std::ostream & operator<<(std::ostream &ost, const MPQCData &data);

src/Fragmentation/Summation/Containers/MPQCDataFused.hpp

@@ -22 +22 @@
  * instance when going throughb the list with boost::mpl::for_each.
  */
-namespace MPQCDataFused {
+struct MPQCDataFused {
   // keys for energy_t
   struct energy_total {};
@@ -51 +51 @@
   struct times_gather_cputime {};
   struct times_gather_flops {};
-}
+};
 
 

src/Fragmentation/Summation/Containers/VMGData.cpp

@@ -44 +44 @@
   sampled_potential(_props),
   nuclei_long(0.),
-  electron_long(0.)
+  electron_long(0.),
+  hasForces(false)
 {}
 
@@ -54 +55 @@
   ost << "Nuclei long-Range energy: " << data.nuclei_long << std::endl;
   ost << "Electron long-Range energy: " << data.electron_long << std::endl;
+  ost << "Nuclei long-Range forces: " << data.forces << std::endl;
   return ost;
 }

src/Fragmentation/Summation/Containers/VMGData.hpp

@@ -22 +22 @@
 #include <vector>
 
+#include "Fragmentation/Summation/SetValues/FragmentForces.hpp"
 #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
 
@@ -48 +49 @@
   //!> electron long-range contribution to energy
   double electron_long;
+  //!> force vectors of all nuclei
+  FragmentForces forces;
+
+  //!> internal value to check when we should sum forces and when not
+  bool hasForces;
 
 private:
@@ -61 +67 @@
       ar & nuclei_long;
     ar & electron_long;
+    if (version > 2) {
+      ar & forces;
+      hasForces = true;
+    } else
+      hasForces = false;
   }
 };
 
-BOOST_CLASS_VERSION(VMGData, 2)
+BOOST_CLASS_VERSION(VMGData, 3)
 
 std::ostream & operator<<(std::ostream &ost, const VMGData &data);

src/Fragmentation/Summation/Containers/VMGDataFused.hpp

@@ -23 +23 @@
  * instance when going throughb the list with boost::mpl::for_each.
  */
-namespace VMGDataFused {
+struct VMGDataFused {
   // keys for sampled_potential
   struct sampled_potential {};
@@ -29 +29 @@
   struct nuclei_long {};
   struct electron_long {};
+  // keys for forces
+  struct forces {};
 
   // keys for longrange
-
   struct electron_longrange {};
   struct electron_shortrange {};
@@ -40 +41 @@
   struct total_longrange {};
   struct total_shortrange {};
-}
+
+  struct forces_longrange {};
+  struct forces_shortrange {};
+};
 
 

src/Fragmentation/Summation/Containers/VMGDataMap.hpp

@@ -18 +18 @@
 #include <boost/mpl/list.hpp>
 
+#include <vector>
+
 #include "Fragmentation/Summation/Containers/VMGDataFused.hpp"
 
 class SamplingGrid;
+class IndexedVectors;
 
 /** This boost::fusion map defines key-value or rather key-type pairs with
@@ -30 +33 @@
  */
 typedef boost::fusion::map<
-    boost::fusion::pair<VMGDataFused::sampled_potential, SamplingGrid >,
-    boost::fusion::pair<VMGDataFused::both_sampled_potential, SamplingGrid >,
     boost::fusion::pair<VMGDataFused::nuclei_long, double >,
     boost::fusion::pair<VMGDataFused::electron_long, double >
@@ -37 +38 @@
 
 typedef boost::mpl::list<
-    VMGDataFused::sampled_potential,
-    VMGDataFused::both_sampled_potential,
     VMGDataFused::nuclei_long,
     VMGDataFused::electron_long
 > VMGDataVector_t;
+
+/** This boost::fusion map defines key-value or rather key-type pairs with
+ * which we associate all forces data members in VMGData and their type.
+ *
+ * This lets us resolves any ambiguitites of types in VMGData, e.g.
+ * to know vector<double> is forces or energy_eigenvalues.
+ *
+ */
+typedef boost::fusion::map<
+    boost::fusion::pair<VMGDataFused::forces, IndexedVectors >
+> VMGDataForceMap_t;
+
+typedef boost::mpl::list<
+    VMGDataFused::forces
+> VMGDataForceVector_t;
+
+typedef boost::fusion::map<
+    boost::fusion::pair<VMGDataFused::sampled_potential, SamplingGrid >,
+    boost::fusion::pair<VMGDataFused::both_sampled_potential, SamplingGrid >
+> VMGDataGridMap_t;
+
+typedef boost::mpl::list<
+    VMGDataFused::sampled_potential,
+    VMGDataFused::both_sampled_potential
+> VMGDataGridVector_t;
 
 /** This boost::fusion map defines key-value or rather key-type pairs with
@@ -72 +96 @@
 > VMGDataLongRangeVector_t;
 
+typedef boost::fusion::map<
+    boost::fusion::pair<VMGDataFused::forces_longrange, IndexedVectors >,
+    boost::fusion::pair<VMGDataFused::forces_shortrange, IndexedVectors >
+> VMGDataLongRangeForceMap_t;
+
+typedef boost::mpl::list<
+    VMGDataFused::forces_longrange,
+    VMGDataFused::forces_shortrange
+> VMGDataLongRangeForceVector_t;
+
 
 #endif /* VMGDATAMAP_HPP_ */

src/Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp

@@ -38 +38 @@
   (nuclei_long) \
   (electron_long) \
+  (forces) \
   (electron_longrange) \
   (electron_shortrange) \
@@ -45 +46 @@
   (nuclei_shortrange) \
   (total_longrange) \
-  (total_shortrange)
+  (total_shortrange) \
+  (forces_longrange) \
+  (forces_shortrange)
 
 /// we take note of the number of keys in tokensequence as (local) loop bounds below

src/Fragmentation/Summation/Converter/DataConverter.hpp

@@ -88 +88 @@
     LOG(4, "DEBUG: Current extracted Data is " << extractedData << ".");
     VMGDataMap_t instance;
-    boost::fusion::at_key<VMGDataFused::sampled_potential>(instance) = extractedData.sampled_potential;
-    boost::fusion::at_key<VMGDataFused::both_sampled_potential>(instance) = extractedData.both_sampled_potential;
     boost::fusion::at_key<VMGDataFused::nuclei_long>(instance) = extractedData.nuclei_long;
     boost::fusion::at_key<VMGDataFused::electron_long>(instance) = extractedData.electron_long;
@@ -95 +93 @@
   }
 }
-#endif
-
-inline void convertMPQCDatatoForceMap(
-    const std::map<JobId_t, MPQCData> &fragmentData,
+
+template <>
+inline void convertDataTo<VMGData, VMGDataGridMap_t>(
+    const std::map<JobId_t, VMGData> &longrangeData,
+    std::map<JobId_t, VMGDataGridMap_t> &VMGData_fused)
+{
+  // energy_t
+  VMGData_fused.clear();
+  for(std::map<JobId_t, VMGData>::const_iterator dataiter = longrangeData.begin();
+      dataiter != longrangeData.end(); ++dataiter) {
+    const VMGData &extractedData = dataiter->second;
+    LOG(4, "DEBUG: Current extracted Data is " << extractedData << ".");
+    VMGDataGridMap_t instance;
+    boost::fusion::at_key<VMGDataFused::sampled_potential>(instance) = extractedData.sampled_potential;
+    boost::fusion::at_key<VMGDataFused::both_sampled_potential>(instance) = extractedData.both_sampled_potential;
+    VMGData_fused.insert( std::make_pair(dataiter->first, instance) );
+  }
+}
+#endif
+
+template <class datatype, class dataforcemap, class dataforcefused>
+inline void convertDatatoForceMap(
+    const std::map<JobId_t, datatype> &fragmentData,
     const KeySetsContainer &ForceKeySet,
-    std::map<JobId_t, MPQCDataForceMap_t> &MPQCData_Force_fused)
+    std::map<JobId_t, dataforcemap> &Data_Force_fused)
 {
   // forces
   ASSERT( ForceKeySet.KeySets.size() == fragmentData.size(),
       "FragmentationAutomationAction::performCall() - indices and fragmentData differ in size.");
-  MPQCData_Force_fused.clear();
-  std::map<JobId_t, MPQCData>::const_iterator dataiter = fragmentData.begin();
+  Data_Force_fused.clear();
+  typename std::map<JobId_t, datatype>::const_iterator dataiter = fragmentData.begin();
   KeySetsContainer::ArrayOfIntVectors::const_iterator arrayiter = ForceKeySet.KeySets.begin();
   for(;dataiter != fragmentData.end(); ++dataiter, ++arrayiter) {
-    const MPQCData &extractedData = dataiter->second;
-    LOG(4, "DEBUG: Current extracted Data is " << extractedData << ".");
-    MPQCDataForceMap_t instance;
+    const datatype &extractedData = dataiter->second;
+    LOG(4, "DEBUG: Current extracted Data is " << extractedData << ".");
+    dataforcemap instance;
     // must convert int to index_t
     if (DoLog(5)) {
@@ -122 +139 @@
     }
     IndexedVectors::indices_t indices(arrayiter->begin(), arrayiter->end());
-    boost::fusion::at_key<MPQCDataFused::forces>(instance) =
+    boost::fusion::at_key<typename dataforcefused::forces>(instance) =
         IndexedVectors(indices, extractedData.forces);
-    MPQCData_Force_fused.insert( std::make_pair(dataiter->first, instance) );
+    Data_Force_fused.insert( std::make_pair(dataiter->first, instance) );
   }
 }

src/Fragmentation/Summation/Makefile.am

@@ -77 +77 @@
 libMolecuilderFragmentationSummation_libincludedir = $(libdir)/MoleCuilder/include
 nodist_libMolecuilderFragmentationSummation_libinclude_HEADERS = $(top_builddir)/libmolecuilder_config.h
-
-## Install the generated pkg-config file (.pc) into the expected location for
-## architecture-dependent package configuration information.  Occasionally,
-## pkg-config files are also used for architecture-independent data packages,
-## in which case the correct install location would be $(datadir)/pkgconfig.
-#pkgconfigdir = $(libdir)/pkgconfig
-#pkgconfig_DATA = $(top_builddir)/MoleCuilder.pc

src/Fragmentation/Summation/SetValues/Makefile.am

@@ -5 +5 @@
         Fragmentation/Summation/SetValues/Eigenvalues.cpp \
         Fragmentation/Summation/SetValues/Fragment.cpp \
+        Fragmentation/Summation/SetValues/FragmentForces.cpp \
         Fragmentation/Summation/SetValues/Histogram.cpp \
         Fragmentation/Summation/SetValues/IndexedVectors.cpp \
@@ -13 +14 @@
         Fragmentation/Summation/SetValues/Eigenvalues.hpp \
         Fragmentation/Summation/SetValues/Fragment.hpp \
+        Fragmentation/Summation/SetValues/FragmentForces.hpp \
         Fragmentation/Summation/SetValues/Histogram.hpp \
         Fragmentation/Summation/SetValues/IndexedVectors.hpp \