Changes in src/Fragmentation/Exporters/SaturatedFragment.cpp [5aaa43:5d5550]

File:

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

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()
-{
-  // 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;
-  for (std::vector<atom *>::const_iterator iter = atoms.begin();
-      iter != atoms.end();
+  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 (atoms_t::const_iterator iter = _atoms.begin();
+      iter != _atoms.end();
       ++iter) {
     atom * const Walker = *iter;
@@ -125 +142 @@
         ++BondRunner) {
       atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
-      // if in 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 (!AddHydrogenReplacementAtom(
-              (*BondRunner),
-              Walker,
-              OtherWalker,
-              World::getInstance().getConfig()->IsAngstroem == 1))
-            exit(1);
+          if (CutBonds.count(Walker) == 0)
+            CutBonds.insert( std::make_pair(Walker, BondList() ));
+          CutBonds[Walker].push_back(*BondRunner);
+      }
+    }
+  }
+
+  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());
+  }
+
+  // go through all cut bonds and replace with a hydrogen
+  for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
+      atomiter != CutBonds.end(); ++atomiter) {
+    atom * const Walker = atomiter->first;
+    if (!saturateAtom(Walker, atomiter->second))
+      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.;
+//          }
+//        }
+
+        // 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, 1., father);
+          FullMolecule.insert(hydrogen.getId());
         }
-//        } 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();
       }
     }
-  }
+  } 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;
+}
+
+bool SaturatedFragment::saturateAtom(
+    atom * const _atom,
+    const BondList &_cutbonds)
+{
+  // go through each bond and replace
+  for (BondList::const_iterator bonditer = _cutbonds.begin();
+      bonditer != _cutbonds.end(); ++bonditer) {
+    atom * const OtherWalker = (*bonditer)->GetOtherAtom(_atom);
+    if (!AddHydrogenReplacementAtom(
+        (*bonditer),
+        _atom,
+        OtherWalker,
+        World::getInstance().getConfig()->IsAngstroem == 1))
+      return false;
+  }
+  return true;
 }
 
@@ -270 +444 @@
   if (BondRescale == -1) {
     ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
-    return false;
-    BondRescale = bondlength;
+    BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree());
+    if (BondRescale == -1) {
+      ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of any degree!");
+      return false;
+      BondRescale = bondlength;
+    }
   } else {
     if (!IsAngstroem)