Changeset f5fa48

Timestamp:

Aug 20, 2014, 1:06:47 PM (11 years ago)

Author:

Frederik Heber <heber@…>

Children:

d635829

Parents:

ce0ca4

git-author:

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

git-committer:

Frederik Heber <heber@…> (08/20/14 13:06:47)

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) (8 diffs)
• Fragmentation/Exporters/SaturatedFragment.hpp (modified) (5 diffs)
• Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp (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

@@ -65 +65 @@
     HydrogenPool &_hydrogens,
     const enum HydrogenTreatment _treatment,
-    const enum HydrogenSaturation _saturation) :
+    const enum HydrogenSaturation _saturation,
+    const GlobalSaturationPositions_t &_globalsaturationpositions) :
   container(_container),
   set(_set),
@@ -79 +80 @@
   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);
 }
 
@@ -101 +106 @@
 }
 
-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);
@@ -117 +122 @@
   LOG(4, "DEBUG: We have gathered the following atoms: " << atoms);
 
-//  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,
+    KeySet &_FullMolecule)
+{
+  //  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;
@@ -135 +149 @@
       atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
       // if other atom is in key set
-      if (set.find(OtherWalker->getId()) != set.end()) {
+      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)
@@ -150 +164 @@
         LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
             << *OtherWalker << ", who is not in this fragment molecule.");
-        if ((treatment == ExcludeHydrogen) && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
+        if ((_treatment == ExcludeHydrogen) && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
           LOG(4, "REJECT: " << *OtherWalker << " is a hydrogen, that are excluded from the set.");
-          FullMolecule.insert(OtherWalker->getId());
+          _FullMolecule.insert(OtherWalker->getId());
         } else {
           LOG(3, "ACCEPT: Adding " << **BondRunner << " as a cut bond.");
@@ -162 +176 @@
   }
   LOG(4, "DEBUG: We have gathered the following CutBonds: " << CutBonds);
+
+  return CutBonds;
+}
+
+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, FullMolecule);
 
   // go through all cut bonds and replace with a hydrogen
@@ -172 +198 @@
       if (!saturateAtom(Walker, atomiter->second))
         exit(1);
+    }
+  } else
+    LOG(3, "INFO: We are not saturating cut bonds.");
+}
+
+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, FullMolecule);
+
+  // 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

src/Fragmentation/Exporters/SaturatedFragment.hpp

@@ -23 +23 @@
 #include "Parser/FormatParserStorage.hpp"
 
+#include "LinearAlgebra/Vector.hpp"
+
 class atom;
 class HydrogenPool;
-class Vector;
 
 /** The SaturatedFragment class acts as a wrapper to a KeySet by adding a list
@@ -43 +44 @@
   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 +57 @@
    * \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 +64 @@
       HydrogenPool &_hydrogens,
       const enum HydrogenTreatment _treatment,
-      const enum HydrogenSaturation saturation);
+      const enum HydrogenSaturation saturation,
+      const GlobalSaturationPositions_t &_globalsaturationpositions);
 
   /** Destructor of class SaturatedFragment.
@@ -101 +111 @@
   /** 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.
+   *
+   * \param _atom atom whose cut bonds to saturate
+   * \param _cutbonds list of cut bonds for \a _atom
+   * \return true - bonds saturated, false - something went wrong
+   */
+  bool saturateAtom(atom * const _atom, const BondList &_cutbonds);
 
   /** Helper function to get a hydrogen replacement for a given \a replacement.

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);
 }