Changeset e139180

Timestamp:

Aug 20, 2014, 1:04:08 PM (11 years ago)

Author:

Frederik Heber <heber@…>

Children:

d734ff

Parents:

0d5ca7

git-author:

Frederik Heber <heber@…> (06/04/14 11:22:49)

git-committer:

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

Message:

Fixes to SaturatedFragment::saturateFragment().

• properly setting up the number of points and the "old" polygon.
• properly filling in the hydrogen atoms at the calculated places.
• We have the number of remaining bonds plus the rest. The rest is the valence minus the number of remaining bonds each weighted with its degree. This gives the right number of places to put hydrogens and fill up the valence.
• TESTS: Removed XFAIL from FragmentMolecule cycles regression test.

Files:

• src/Fragmentation/Exporters/SaturatedFragment.cpp (modified) (8 diffs)
• src/Fragmentation/Exporters/SaturatedFragment.hpp (modified) (1 diff)
• tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at (modified) (1 diff)

src/Fragmentation/Exporters/SaturatedFragment.cpp ¶

@@ -38 +38 @@
 #include "SaturatedFragment.hpp"
 
+#include <algorithm>
 #include <cmath>
 #include <iostream>
@@ -114 +115 @@
     atoms.push_back(Walker);
   }
+  LOG(4, "DEBUG: We have gathered the following atoms: " << atoms);
 
 //  bool LonelyFlag = false;
@@ -123 +125 @@
       ++iter) {
     atom * const Walker = *iter;
+    // start with an empty list
+    CutBonds.insert( std::make_pair(Walker, BondList() ));
 
     // go through all bonds
@@ -146 +150 @@
         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() ));
+        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());
+        } else {
+          LOG(3, "ACCEPT: Adding " << **BondRunner << " as a cut bond.");
+          // there is always at least an empty list
           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();
       }
     }
   }
+  LOG(4, "DEBUG: We have gathered the following CutBonds: " << CutBonds);
 
   // 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);
-  }
+  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 " << *Walker);
+
+      if (!saturateAtom(Walker, atomiter->second))
+        exit(1);
+    }
+  } else
+    LOG(3, "INFO: We are not saturating cut bonds.");
 }
 
@@ -188 +194 @@
 //  }
 
-  // gather the nodes of the shape defined by the current set of bonded atoms
   SphericalPointDistribution::Polygon_t Polygon;
-  for (BondList::const_iterator bonditer = _cutbonds.begin();
-      bonditer != _cutbonds.end(); ++bonditer) {
-    Vector DistanceVector;
-    if ((*bonditer)->leftatom == _atom)
-      DistanceVector = (*bonditer)->rightatom->getPosition() - (*bonditer)->leftatom->getPosition();
-    else
-      DistanceVector = (*bonditer)->leftatom->getPosition() - (*bonditer)->rightatom->getPosition();
-    // always use unit distances
-    DistanceVector.Normalize();
-    Polygon.push_back( DistanceVector );
-  }
-  LOG(3, "DEBUG: Polygon of atom " << _atom << " to saturate is " << Polygon);
-
-  // get the new number of bonds (where all cut bonds are replaced by as
-  // many hydrogens as is their degree)
-  //  unsigned int NumberOfPoints = _atom->getElement().NoValenceOrbitals;
-  unsigned int NumberOfPoints =
-      _atom->getListOfBonds().size() - _cutbonds.size();
-  const BondList& ListOfBonds = _atom->getListOfBonds();
-  for (BondList::const_iterator BondRunner = ListOfBonds.begin();
-      BondRunner != ListOfBonds.end();
-      ++BondRunner) {
-    NumberOfPoints += (*BondRunner)->getDegree();
-    // ALTERNATIVE starting from valence:
-    //    // if not one of the cut bonds, reduce by its bond degree -1 (for the one bond itself)
-    //    if (_cutbonds.count(*BondRunner))
-    //      NumberOfPoints -= (*BondRunner)->getDegree() - 1;
-  }
-  LOG(3, "DEBUG: There are " << NumberOfPoints
-      << " to fill in in total for this atom " << _atom << ".");
+  {
+    // prepare a list of "uncut" bonds via set_difference. For this both lists
+    // have to be sorted.
+    typedef std::vector<bond::ptr> BondVector_t;
+    BondVector_t ListOfBonds(_atom->getListOfBonds().begin(),_atom->getListOfBonds().end());
+    std::sort(ListOfBonds.begin(), ListOfBonds.end());
+    BondVector_t CutBonds(_cutbonds.begin(), _cutbonds.end());
+    std::sort(CutBonds.begin(), CutBonds.end());
+    const BondVector_t::iterator eraseiter = std::set_difference(
+        ListOfBonds.begin(), ListOfBonds.end(),
+        CutBonds.begin(), CutBonds.end(),
+        ListOfBonds.begin());
+    ListOfBonds.erase(eraseiter, ListOfBonds.end());
+
+    // gather the nodes of the shape defined by the current set of bonded atoms
+    for (BondVector_t::const_iterator bonditer = ListOfBonds.begin();
+        bonditer != ListOfBonds.end();
+        ++bonditer) {
+      Vector DistanceVector;
+      if ((*bonditer)->leftatom == _atom)
+        DistanceVector = (*bonditer)->rightatom->getPosition() - (*bonditer)->leftatom->getPosition();
+      else
+        DistanceVector = (*bonditer)->leftatom->getPosition() - (*bonditer)->rightatom->getPosition();
+      // always use unit distances
+      DistanceVector.Normalize();
+      Polygon.push_back( DistanceVector );
+    }
+    LOG(3, "DEBUG: Polygon of atom " << *_atom << " to saturate is " << Polygon);
+  }
+
+  unsigned int NumberOfPoints = 0;
+  {
+    // get the new number of bonds
+    const BondList & ListOfBonds = _atom->getListOfBonds();
+    NumberOfPoints = ListOfBonds.size() - _cutbonds.size(); // number remaining bonds
+    NumberOfPoints += _atom->getElement().getNoValenceOrbitals(); // plus valence
+    // minus remaining bonds weighted by its degree
+    for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner) {
+      // if not one of the cut bonds, reduce by its bond degree -1 (for the one bond itself)
+      const BondList::const_iterator finditer =
+          std::find(_cutbonds.begin(), _cutbonds.end(), *BondRunner);
+      if (finditer == _cutbonds.end())
+        NumberOfPoints -= (*BondRunner)->getDegree();
+    }
+    LOG(3, "DEBUG: There are " << NumberOfPoints
+        << " places to fill in in total for this atom " << *_atom << ".");
+  }
 
   // get perfect node distribution for the given remaining atoms with respect
@@ -287 +311 @@
     atom * const _father = _atom;
     LOG(4, "DEBUG: Filling saturation hydrogen for atom " << _atom << " at " << *iter);
-    setHydrogenReplacement(
+    const atom& hydrogen = setHydrogenReplacement(
         _atom,
         *iter,
         1.,
         _father);
+    FullMolecule.insert(hydrogen.getId());
   }
 
@@ -355 +380 @@
 }
 
-void SaturatedFragment::setHydrogenReplacement(
+const atom& SaturatedFragment::setHydrogenReplacement(
     const atom * const _OwnerAtom,
     const Vector &_position,
@@ -368 +393 @@
   _atom->father = _father;
   SaturationHydrogens.insert(_atom->getId());
+  return *_atom;
 }
 

src/Fragmentation/Exporters/SaturatedFragment.hpp ¶

@@ -117 +117 @@
    * \param _distance scale factor to the distance (default 1.)
    * \param _father bond partner of \a _OwnerAtom that is replaced
+   * \return pointer to saturation hydrogen atom
    */
-  void setHydrogenReplacement(
+  const atom& setHydrogenReplacement(
       const atom * const _OwnerAtom,
       const Vector &_position,

tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at ¶

@@ -37 +37 @@
 AT_SETUP([Fragmentation - Fragmentation with cycles])
 AT_KEYWORDS([fragmentation fragment-molecule cycle])
-AT_XFAIL_IF([/bin/true])
 
 file=benzene.pdb