Changeset 98a293b

Timestamp:

Aug 16, 2014, 3:52:30 PM (10 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:

a8f1d2

Parents:

bdd5a8

git-author:

Frederik Heber <heber@…> (07/18/14 17:08:09)

git-committer:

Frederik Heber <heber@…> (08/16/14 15:52:30)

Message:

SaturatedFragment can deal with a global saturation position map.

• so far, we create an empty one in FragmentationAction such that nothing's changed for the moment.
• similarly in StoreSaturatedFragmentAction. However, there this is intended as only local information is required (it's only a single fragment).

Location:

src

Files:

• Actions/FragmentationAction/FragmentationAction.cpp (modified) (4 diffs)
• Actions/FragmentationAction/StoreSaturatedFragmentAction.cpp (modified) (2 diffs)
• Fragmentation/Exporters/ExportGraph.cpp (modified) (2 diffs)
• Fragmentation/Exporters/ExportGraph.hpp (modified) (2 diffs)
• Fragmentation/Exporters/ExportGraph_ToFiles.cpp (modified) (1 diff)
• Fragmentation/Exporters/ExportGraph_ToFiles.hpp (modified) (1 diff)
• Fragmentation/Exporters/ExportGraph_ToJobs.cpp (modified) (1 diff)
• Fragmentation/Exporters/ExportGraph_ToJobs.hpp (modified) (1 diff)
• Fragmentation/Exporters/SaturatedFragment.cpp (modified) (5 diffs)
• Fragmentation/Exporters/SaturatedFragment.hpp (modified) (6 diffs)
• Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp (modified) (1 diff)
• Fragmentation/Summation/Converter/DataConverter.hpp (modified) (1 diff)

src/Actions/FragmentationAction/FragmentationAction.cpp

@@ -42 +42 @@
 #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
 #include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
+#include "Fragmentation/Exporters/SaturatedFragment.hpp"
 #include "Fragmentation/Fragmentation.hpp"
 #include "Fragmentation/Graph.hpp"
@@ -261 +262 @@
   }
 
+  // create global saturation positions map
+  SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
+
   {
     const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
@@ -266 +270 @@
     if (params.types.get().size() != 0) {
       // store molecule's fragment to file
-      ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation);
+      ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
       exporter.setPrefix(params.prefix.get());
       exporter.setOutputTypes(params.types.get());
@@ -272 +276 @@
     } else {
       // store molecule's fragment in FragmentJobQueue
-      ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation);
+      ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
       exporter.setLevel(params.level.get());
       exporter();

src/Actions/FragmentationAction/StoreSaturatedFragmentAction.cpp

@@ -39 +39 @@
 #include "CodePatterns/Log.hpp"
 #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
+#include "Fragmentation/Exporters/SaturatedFragment.hpp"
 #include "Fragmentation/Graph.hpp"
 #include "World.hpp"
@@ -77 +78 @@
   // store molecule's fragment to file
   {
+    // we use an empty map here such that saturation is done locally
+    SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
+
     const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
-    ExportGraph_ToFiles exporter(TotalGraph, IncludeHydrogen, saturation);
+    ExportGraph_ToFiles exporter(TotalGraph, IncludeHydrogen, saturation, globalsaturationpositions);
     exporter.setPrefix(params.prefix.get());
     exporter.setOutputTypes(params.types.get());

src/Fragmentation/Exporters/ExportGraph.cpp

@@ -58 +58 @@
   const Graph &_graph,
   const enum HydrogenTreatment _treatment,
-  const enum HydrogenSaturation _saturation) :
+  const enum HydrogenSaturation _saturation,
+  const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
   TotalGraph(_graph),
   BondFragments(World::getPointer()),
   treatment(_treatment),
   saturation(_saturation),
+  globalsaturationpositions(_globalsaturationpositions),
   CurrentKeySet(TotalGraph.begin())
 {
@@ -115 +117 @@
           hydrogens,
           treatment,
-          saturation)
+          saturation,
+          globalsaturationpositions)
   );
   // and return

src/Fragmentation/Exporters/ExportGraph.hpp

@@ -41 +41 @@
       const Graph &_graph,
       const enum HydrogenTreatment _treatment,
-      const enum HydrogenSaturation _saturation);
+      const enum HydrogenSaturation _saturation,
+      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
   virtual ~ExportGraph();
 
@@ -102 +103 @@
   const enum HydrogenSaturation saturation;
 
+  //!> Global information over all atoms with saturation positions to be used per fragment
+  const SaturatedFragment::GlobalSaturationPositions_t &globalsaturationpositions;
+
 private:
   //!> iterator pointing at the CurrentKeySet to be exported

src/Fragmentation/Exporters/ExportGraph_ToFiles.cpp

@@ -57 +57 @@
  * @param _treatment whether to always add already present hydrogens or not
  * @param _saturation whether to saturate dangling bonds with hydrogen or not
+ * @param _globalsaturationpositions possibly empty map with global information
+ *        where to place saturation hydrogens to fulfill consistency principle
  */
 ExportGraph_ToFiles::ExportGraph_ToFiles(
     const Graph &_graph,
     const enum HydrogenTreatment _treatment,
-    const enum HydrogenSaturation _saturation) :
-                ExportGraph(_graph, _treatment, _saturation)
+    const enum HydrogenSaturation _saturation,
+    const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
+                ExportGraph(_graph, _treatment, _saturation, _globalsaturationpositions)
 {}
 

src/Fragmentation/Exporters/ExportGraph_ToFiles.hpp

@@ -33 +33 @@
             const Graph &_graph,
             const enum HydrogenTreatment _treatment,
-            const enum HydrogenSaturation _saturation);
+            const enum HydrogenSaturation _saturation,
+      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
         virtual ~ExportGraph_ToFiles();
 

src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

@@ -59 +59 @@
     const Graph &_graph,
     const enum HydrogenTreatment _treatment,
-    const enum HydrogenSaturation _saturation) :
-    ExportGraph(_graph, _treatment, _saturation),
+    const enum HydrogenSaturation _saturation,
+    const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
+    ExportGraph(_graph, _treatment, _saturation,_globalsaturationpositions),
     level(5)
 {}

src/Fragmentation/Exporters/ExportGraph_ToJobs.hpp

@@ -35 +35 @@
    * \param _treatment whether hydrogen is excluded in the _graph or not
    * \param _saturation whether we saturate dangling bonds or not
+   * \param _globalsaturationpositions possibly empty map with global information
+   *        where to place saturation hydrogens to fulfill consistency principle
    */
   ExportGraph_ToJobs(
       const Graph &_graph,
       const enum HydrogenTreatment _treatment,
-      const enum HydrogenSaturation _saturation);
+      const enum HydrogenSaturation _saturation,
+      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
   virtual ~ExportGraph_ToJobs();
 

src/Fragmentation/Exporters/SaturatedFragment.cpp

@@ -63 +63 @@
     HydrogenPool &_hydrogens,
     const enum HydrogenTreatment _treatment,
-    const enum HydrogenSaturation _saturation) :
+    const enum HydrogenSaturation _saturation,
+    const GlobalSaturationPositions_t &_globalsaturationpositions) :
   container(_container),
   set(_set),
@@ -77 +78 @@
   container.insert(set);
 
-  // prepare saturation hydrogens
-  saturate();
+  // prepare saturation hydrogens, either using global information
+  // or if not given, local information (created in the function)
+  if (_globalsaturationpositions.empty())
+    saturate();
+  else
+    saturate(_globalsaturationpositions);
 }
 
@@ -99 +104 @@
 }
 
-void SaturatedFragment::saturate()
-{
-  // so far, we just have a set of keys. Hence, convert to atom refs
-  // and gather all atoms in a vector
-  std::vector<atom *> atoms;
-  for (KeySet::const_iterator iter = FullMolecule.begin();
-      iter != FullMolecule.end();
+typedef std::vector<atom *> atoms_t;
+
+atoms_t gatherAllAtoms(const KeySet &_FullMolecule)
+{
+  atoms_t atoms;
+  for (KeySet::const_iterator iter = _FullMolecule.begin();
+      iter != _FullMolecule.end();
       ++iter) {
     atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
     ASSERT( Walker != NULL,
-        "SaturatedFragment::OutputConfig() - id "
+        "gatherAllAtoms() - id "
         +toString(*iter)+" is unknown to World.");
     atoms.push_back(Walker);
   }
 
-//  bool LonelyFlag = false;
-  // go through each atom of the fragment and gather all cut bonds in list
-  typedef std::map<atom *, BondList > CutBonds_t;
+  return atoms;
+}
+
+typedef std::map<atom *, BondList > CutBonds_t;
+
+CutBonds_t gatherCutBonds(
+    const atoms_t &_atoms,
+    const KeySet &_set,
+    const enum HydrogenTreatment _treatment)
+{
+  //  bool LonelyFlag = false;
   CutBonds_t CutBonds;
-  for (std::vector<atom *>::const_iterator iter = atoms.begin();
-      iter != atoms.end();
+  for (atoms_t::const_iterator iter = _atoms.begin();
+      iter != _atoms.end();
       ++iter) {
     atom * const Walker = *iter;
@@ -129 +142 @@
         ++BondRunner) {
       atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
-      // if other atom is in key set
-      if (set.find(OtherWalker->getId()) != set.end()) {
+      // if other atom is in key set or is a specially treated hydrogen
+      if (_set.find(OtherWalker->getId()) != _set.end()) {
         LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *OtherWalker << ".");
-//        if (OtherWalker->getId() > Walker->getId()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
-////          std::stringstream output;
-////            output << "ACCEPT: Adding Bond: "
-//          output << Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
-////            LOG(3, output.str());
-//          //NumBonds[(*iter)->getNr()]++;
-//        } else {
-////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
-//        }
-//        LonelyFlag = false;
+      } else if ((_treatment == ExcludeHydrogen)
+          && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
+        LOG(4, "DEBUG: Walker " << *Walker << " is bound to specially treated hydrogen " <<
+            *OtherWalker << ".");
       } else {
         LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
             << *OtherWalker << ", who is not in this fragment molecule.");
-        if (saturation == DoSaturate) {
-//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
           if (CutBonds.count(Walker) == 0)
             CutBonds.insert( std::make_pair(Walker, BondList() ));
           CutBonds[Walker].push_back(*BondRunner);
-        }
-//        } else if ((treatment == ExcludeHydrogen) && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
-//          // just copy the atom if it's a hydrogen
-//          atom * const OtherWalker = Leaf->AddCopyAtom(OtherWalker);
-//          Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
-//        }
-        //NumBonds[(*iter)->getNr()] += Binder->getDegree();
       }
     }
+  }
+
+  return CutBonds;
+}
+
+typedef std::vector<atomId_t> atomids_t;
+
+atomids_t gatherPresentExcludedHydrogens(
+    const atoms_t &_atoms,
+    const KeySet &_set,
+    const enum HydrogenTreatment _treatment)
+{
+  //  bool LonelyFlag = false;
+  atomids_t ExcludedHydrogens;
+  for (atoms_t::const_iterator iter = _atoms.begin();
+      iter != _atoms.end();
+      ++iter) {
+    atom * const Walker = *iter;
+
+    // go through all bonds
+    const BondList& ListOfBonds = Walker->getListOfBonds();
+    for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner) {
+      atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
+      // if other atom is in key set or is a specially treated hydrogen
+      if (_set.find(OtherWalker->getId()) != _set.end()) {
+        LOG(6, "DEBUG: OtherWalker " << *OtherWalker << " is in set already.");
+      } else if ((_treatment == ExcludeHydrogen)
+          && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
+        LOG(5, "DEBUG: Adding excluded hydrogen OtherWalker " << *OtherWalker << ".");
+        ExcludedHydrogens.push_back(OtherWalker->getId());
+      } else {
+        LOG(6, "DEBUG: OtherWalker " << *Walker << " is not in this fragment molecule and no hydrogen.");
+      }
+    }
+  }
+
+  return ExcludedHydrogens;
+}
+
+void SaturatedFragment::saturate()
+{
+  // so far, we just have a set of keys. Hence, convert to atom refs
+  // and gather all atoms in a vector
+  std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
+
+  // go through each atom of the fragment and gather all cut bonds in list
+  CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
+
+  // add excluded hydrogens to FullMolecule if treated specially
+  if (treatment == ExcludeHydrogen) {
+    atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
+    FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
   }
 
@@ -168 +220 @@
       exit(1);
   }
+}
+
+void SaturatedFragment::saturate(
+    const GlobalSaturationPositions_t &_globalsaturationpositions)
+{
+  // so far, we just have a set of keys. Hence, convert to atom refs
+  // and gather all atoms in a vector
+  std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
+
+  // go through each atom of the fragment and gather all cut bonds in list
+  CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
+
+  // add excluded hydrogens to FullMolecule if treated specially
+  if (treatment == ExcludeHydrogen) {
+    atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
+    FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
+  }
+
+  // go through all cut bonds and replace with a hydrogen
+  if (saturation == DoSaturate) {
+    for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
+        atomiter != CutBonds.end(); ++atomiter) {
+      atom * const Walker = atomiter->first;
+      LOG(4, "DEBUG: We are now saturating dangling bonds of " << *Walker);
+
+      // gather set of positions for this atom from global map
+      GlobalSaturationPositions_t::const_iterator mapiter =
+          _globalsaturationpositions.find(Walker->getId());
+      ASSERT( mapiter != _globalsaturationpositions.end(),
+          "SaturatedFragment::saturate() - no global information for "
+          +toString(*Walker));
+      const SaturationsPositionsPerNeighbor_t &saturationpositions =
+          mapiter->second;
+
+      // go through all cut bonds for this atom
+      for (BondList::const_iterator bonditer = atomiter->second.begin();
+          bonditer != atomiter->second.end(); ++bonditer) {
+        atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
+
+        // get positions from global map
+        SaturationsPositionsPerNeighbor_t::const_iterator positionsiter =
+            saturationpositions.find(OtherWalker->getId());
+        ASSERT(positionsiter != saturationpositions.end(),
+            "SaturatedFragment::saturate() - no information on bond neighbor "
+            +toString(*OtherWalker)+" to atom "+toString(*Walker));
+        ASSERT(!positionsiter->second.empty(),
+            "SaturatedFragment::saturate() - no positions for saturating bond to"
+            +toString(*OtherWalker)+" to atom "+toString(*Walker));
+
+        // get typical bond distance from elements database
+        double BondDistance = Walker->getType()->getHBondDistance(positionsiter->second.size()-1);
+        if (BondDistance < 0.) {
+          ELOG(2, "saturateAtoms() - no typical hydrogen bond distance of degree "
+              +toString(positionsiter->second.size())+" for element "
+              +toString(Walker->getType()->getName()));
+          // try bond degree 1 distance
+          BondDistance = Walker->getType()->getHBondDistance(1-1);
+          if (BondDistance < 0.) {
+            ELOG(1, "saturateAtoms() - no typical hydrogen bond distance for element "
+                +toString(Walker->getType()->getName()));
+            BondDistance = 1.;
+          }
+        }
+        ASSERT( BondDistance > 0.,
+            "SaturatedFragment::saturate() - negative bond distance");
+
+        // place hydrogen at each point
+        LOG(4, "DEBUG: Places to saturate for atom " << *OtherWalker
+            << " are " << positionsiter->second);
+        atom * const father = Walker;
+        for (SaturationsPositions_t::const_iterator positer = positionsiter->second.begin();
+            positer != positionsiter->second.end(); ++positer) {
+          const atom& hydrogen =
+              setHydrogenReplacement(Walker, *positer, BondDistance, father);
+          FullMolecule.insert(hydrogen.getId());
+        }
+      }
+    }
+  } else
+    LOG(3, "INFO: We are not saturating cut bonds.");
+}
+
+const atom& SaturatedFragment::setHydrogenReplacement(
+  const atom * const _OwnerAtom,
+  const Vector &_position,
+  const double _distance,
+  atom * const _father)
+{
+  atom * const _atom = hydrogens.leaseHydrogen();    // new atom
+  _atom->setPosition( _OwnerAtom->getPosition() + _distance * _position );
+  // always set as fixed ion (not moving during molecular dynamics simulation)
+  _atom->setFixedIon(true);
+  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
+  _atom->father = _father;
+  SaturationHydrogens.insert(_atom->getId());
+
+  return *_atom;
 }
 

src/Fragmentation/Exporters/SaturatedFragment.hpp

@@ -23 +23 @@
 #include "Fragmentation/HydrogenSaturation_enum.hpp"
 #include "Parser/FormatParserStorage.hpp"
+
+#include "LinearAlgebra/Vector.hpp"
 
 class atom;
@@ -43 +45 @@
   typedef std::set<KeySet> KeySetsInUse_t;
 
+  //!> List of points giving saturation positions for a single bond neighbor
+  typedef std::list<Vector> SaturationsPositions_t;
+  //!> map for one atom, containing the saturation points for all its neighbors
+  typedef std::map<int, SaturationsPositions_t> SaturationsPositionsPerNeighbor_t;
+  //!> containing the saturation points over all desired atoms required
+  typedef std::map<int, SaturationsPositionsPerNeighbor_t> GlobalSaturationPositions_t;
+
   /** Constructor of SaturatedFragment requires \a set which we are tightly
    * associated.
@@ -49 +58 @@
    * \param _container container to add KeySet as in-use
    * \param _hydrogens pool with hydrogens for saturation
+   * \param _globalsaturationpositions saturation positions to be used
    */
   SaturatedFragment(
@@ -55 +65 @@
       HydrogenPool &_hydrogens,
       const enum HydrogenTreatment _treatment,
-      const enum HydrogenSaturation saturation);
+      const enum HydrogenSaturation saturation,
+      const GlobalSaturationPositions_t &_globalsaturationpositions);
 
   /** Destructor of class SaturatedFragment.
@@ -101 +112 @@
   /** Helper function to lease and bring in place saturation hydrogens.
    *
+   * Here, we use local information to calculate saturation positions.
+   *
    */
   void saturate();
+
+  /** Helper function to lease and bring in place saturation hydrogens.
+   *
+   * Here, saturation positions have to be calculated before and are fully
+   * stored in \a _globalsaturationpositions.
+   *
+   * \param_globalsaturationpositions
+   */
+  void saturate(const GlobalSaturationPositions_t &_globalsaturationpositions);
 
   /** Replaces all cut bonds with respect to the given atom by hydrogens.
@@ -118 +140 @@
    */
   atom * const getHydrogenReplacement(atom * const replacement);
+
+  /** Sets a saturation hydrogen at the given position in place of \a _father.
+   *
+   * \param _OwnerAtom atom "owning" the hydrogen (i.e. the one getting saturated)
+   * \param _position new position relative to \a _OwnerAtom
+   * \param _distance scale factor to the distance (default 1.)
+   * \param _father bond partner of \a _OwnerAtom that is replaced
+   */
+  const atom& setHydrogenReplacement(
+      const atom * const _OwnerAtom,
+      const Vector &_position,
+      const double _distance,
+      atom * const _father);
 
   /** Leases and adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.

src/Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp

@@ -69 +69 @@
   ASSERT_DO(Assert::Throw);
 
+  SaturatedFragment::GlobalSaturationPositions_t globalpositions;
   set = new KeySet;
-  fragment = new SaturatedFragment(*set, KeySetsInUse, hydrogens, ExcludeHydrogen, DoSaturate);
+  fragment = new SaturatedFragment(
+      *set,
+      KeySetsInUse,
+      hydrogens,
+      ExcludeHydrogen,
+      DoSaturate,
+      globalpositions);
 }
 

src/Fragmentation/Summation/Converter/DataConverter.hpp

@@ -113 +113 @@
     MPQCDataForceMap_t instance;
     // must convert int to index_t
+    if (DoLog(5)) {
+      std::stringstream output;
+      for (KeySetsContainer::IntVector::const_iterator outiter = arrayiter->begin();
+          outiter != arrayiter->end(); ++outiter) {
+        output << *outiter << "\t";
+      }
+      LOG(5, "DEBUG: indices are " << output.str());
+    }
     IndexedVectors::indices_t indices(arrayiter->begin(), arrayiter->end());
     boost::fusion::at_key<MPQCDataFused::forces>(instance) =