Changeset f5fa48 for src/Fragmentation/Exporters/SaturatedFragment.cpp

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).

File:

• src/Fragmentation/Exporters/SaturatedFragment.cpp (modified) (8 diffs)

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