Changeset 97dff0 for src/Fragmentation

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

11cc05

Parents:

972412

git-author:

Frederik Heber <heber@…> (05/29/14 11:29:27)

git-committer:

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

Message:

Reimplemented SaturatedFragment::saturateAtom() and added ::setHydrogenReplacement().

• this is now based on an ideal distribution of electrons around an atom via the VSPER model which we match as closely as possible with the still bound neighboring atoms and fill all remaining vacant spots with hydrogens.
• SaturatedFragment sets father simply to the _given_ atom.
• TEST: dimethyl_bromomalonate marked as XFAIL.

Location:

src/Fragmentation/Exporters

Files:

• SaturatedFragment.cpp (modified) (3 diffs)
• SaturatedFragment.hpp (modified) (2 diffs)

src/Fragmentation/Exporters/SaturatedFragment.cpp

@@ -54 +54 @@
 #include "Descriptors/AtomIdDescriptor.hpp"
 #include "Fragmentation/Exporters/HydrogenPool.hpp"
+#include "Fragmentation/Exporters/SphericalPointDistribution.hpp"
 #include "Fragmentation/HydrogenSaturation_enum.hpp"
 #include "Graph/BondGraph.hpp"
@@ -165 +166 @@
       atomiter != CutBonds.end(); ++atomiter) {
     atom * const Walker = atomiter->first;
-    // go through each bond and replace
-    for (BondList::const_iterator bonditer = atomiter->second.begin();
-        bonditer != atomiter->second.end(); ++bonditer) {
-      atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
-      if (!AddHydrogenReplacementAtom(
-          (*bonditer),
-          Walker,
-          OtherWalker,
-          World::getInstance().getConfig()->IsAngstroem == 1))
-        exit(1);
-    }
-  }
-}
+    if (!saturateAtom(Walker, atomiter->second))
+      exit(1);
+  }
+}
+
+bool SaturatedFragment::saturateAtom(
+    atom * const _atom,
+    const BondList &_cutbonds)
+{
+  // OLD WAY: use AddHydrogenReplacementAtom() on each cut bond
+//  // 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;
+//  }
+
+  // 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 );
+  }
+
+  // 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;
+  }
+
+  // get perfect node distribution for the given remaining atoms with respect
+  // to valence of the atoms (for a saturated fragment, resembles number of bonds)
+  SphericalPointDistribution polygonizer;
+  SphericalPointDistribution::Polygon_t NewPolygon;
+  switch (NumberOfPoints)
+  {
+    case 0:
+      NewPolygon = polygonizer.get<0>();
+      break;
+    case 1:
+      NewPolygon = polygonizer.get<1>();
+      break;
+    case 2:
+      NewPolygon = polygonizer.get<2>();
+      break;
+    case 3:
+      NewPolygon = polygonizer.get<3>();
+      break;
+    case 4:
+      NewPolygon = polygonizer.get<4>();
+      break;
+    case 5:
+      NewPolygon = polygonizer.get<5>();
+      break;
+    case 6:
+      NewPolygon = polygonizer.get<6>();
+      break;
+    case 7:
+      NewPolygon = polygonizer.get<7>();
+      break;
+    case 8:
+      NewPolygon = polygonizer.get<8>();
+      break;
+    case 9:
+      NewPolygon = polygonizer.get<9>();
+      break;
+    case 10:
+      NewPolygon = polygonizer.get<10>();
+      break;
+    case 11:
+      NewPolygon = polygonizer.get<11>();
+      break;
+    case 12:
+      NewPolygon = polygonizer.get<12>();
+      break;
+    case 14:
+      NewPolygon = polygonizer.get<14>();
+      break;
+    default:
+      ASSERT(0, "SaturatedFragment::saturateAtom() - cannot deal with the case "
+          +toString(NumberOfPoints)+".");
+  }
+
+  // then we need to match the old with the new
+  SphericalPointDistribution::Polygon_t RemainingPoints =
+      SphericalPointDistribution::matchSphericalPointDistributions(Polygon, NewPolygon);
+
+  // and place hydrogen atoms at each vacant spot in the distance given by the table
+  for(SphericalPointDistribution::Polygon_t::const_iterator iter = RemainingPoints.begin();
+      iter != RemainingPoints.end(); ++iter) {
+    // find nearest atom as father to this point
+    atom * const _father = _atom;
+    setHydrogenReplacement(
+        _atom,
+        *iter,
+        1.,
+        _father);
+  }
+
+  return true;
+}
+
 
 bool SaturatedFragment::OutputConfig(
@@ -234 +346 @@
   SaturationHydrogens.insert(_atom->getId());
   return _atom;
+}
+
+void 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());
 }
 

src/Fragmentation/Exporters/SaturatedFragment.hpp

@@ -25 +25 @@
 class atom;
 class HydrogenPool;
+class Vector;
 
 /** The SaturatedFragment class acts as a wrapper to a KeySet by adding a list
@@ -110 +111 @@
   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
+   */
+  void 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.
    * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand