Changeset d635829 for src/Actions/FragmentationAction/FragmentationAction.cpp

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

4a44ed, 9d51a8

Parents:

f5fa48

git-author:

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

git-committer:

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

Message:

FragmentationAction now compiles global saturation positions information.

File:

• src/Actions/FragmentationAction/FragmentationAction.cpp (modified) (2 diffs)

src/Actions/FragmentationAction/FragmentationAction.cpp

@@ -43 +43 @@
 #include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
 #include "Fragmentation/Exporters/SaturatedFragment.hpp"
+#include "Fragmentation/Exporters/SphericalPointDistribution.hpp"
 #include "Fragmentation/Fragmentation.hpp"
 #include "Fragmentation/Graph.hpp"
@@ -264 +265 @@
   // create global saturation positions map
   SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
+  {
+    // go through each atom
+    for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
+        iter != world.endAtomSelection(); ++iter) {
+      const atom * const _atom = iter->second;
+
+      // skip hydrogens if treated special
+      const enum HydrogenTreatment treatment =  params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
+      if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {
+        LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);
+        continue;
+      }
+
+      // gather the polygon
+      SphericalPointDistribution::WeightedPolygon_t polygon;
+      const BondList &bondlist = _atom->getListOfBonds();
+      for (BondList::const_iterator bonditer = bondlist.begin();
+          bonditer != bondlist.end(); ++bonditer) {
+        const atom * const _otheratom = (*bonditer)->GetOtherAtom(_atom);
+        Vector bondvector = _otheratom->getPosition() - _atom->getPosition();
+        bondvector.Normalize();
+        polygon.push_back(
+            std::make_pair(
+                bondvector,
+                (*bonditer)->getDegree()
+            ) );
+      }
+      LOG(3, "DEBUG: Polygon is " << polygon << " for this atom " << *_atom << ".");
+
+      // get the valence
+      unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
+      LOG(3, "DEBUG: There are " << NumberOfPoints
+          << " places to fill in in total for this atom " << *_atom << ".");
+
+      // check whether there are any bonds with degree larger than 1
+      unsigned int SumOfDegrees = 0;
+      bool PresentHigherBonds = false;
+      for (BondList::const_iterator bonditer = bondlist.begin();
+          bonditer != bondlist.end(); ++bonditer) {
+        SumOfDegrees += (*bonditer)->getDegree();
+        PresentHigherBonds |= (*bonditer)->getDegree() > 1;
+      }
+
+      // get the association for each bond neighbor
+      SphericalPointDistribution SPD;
+      SphericalPointDistribution::PolygonWithIndexTuples globalinfo;
+      if ((NumberOfPoints != SumOfDegrees) || (PresentHigherBonds))
+        globalinfo = SPD.getAssociatedPoints(polygon, NumberOfPoints);
+      else
+        globalinfo = SphericalPointDistribution::getIdentityAssociation(polygon);
+      LOG(3, "DEBUG: Original polygon was " << polygon << " for this atom " << *_atom << ".");
+      LOG(3, "DEBUG: Got global saturation polygon " << globalinfo.polygon
+          << " for atom " << *_atom);
+      LOG(3, "DEBUG: Associated matching is " << globalinfo.indices
+          << " for atom " << *_atom);
+
+      // convert into the desired entry in the map
+      SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;
+      BondList::const_iterator bonditer = bondlist.begin();
+      // same number of bonds and tuples
+      ASSERT( bondlist.size() == globalinfo.indices.size(),
+          "FragmentationAction::performCall() - number of bonds "
+          +toString(bondlist.size())+" unequal to number of tuples "
+          +toString(globalinfo.indices.size()));
+      for (SphericalPointDistribution::IndexTupleList_t::const_iterator tupleiter =
+          globalinfo.indices.begin();
+          tupleiter != globalinfo.indices.end(); ++tupleiter, ++bonditer) {
+        // note that we cannot check whether the designated vector matches with the
+        // normalized bond vector where it should have originated because the
+        // vector has been rotated.
+//        Vector bondvector =
+//            (*bonditer)->GetOtherAtom(_atom)->getPosition()-_atom->getPosition();
+//        bondvector.Normalize();
+        // no go through all indices of the tuple (and positions)
+        for (SphericalPointDistribution::IndexList_t::const_iterator indexiter =
+            tupleiter->begin(); indexiter != tupleiter->end(); ++indexiter) {
+          // and add them, first trying with new list
+          std::pair<
+              SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,
+              bool
+              > inserter =
+                  positions_per_neighbor.insert(
+                      std::make_pair(
+                          (*bonditer)->GetOtherAtom(_atom)->getId(),
+                          SaturatedFragment::SaturationsPositions_t(
+                              1,
+                              globalinfo.polygon[*indexiter])
+                      )
+                  );
+          // if already pressent, add to this present list
+          if (!inserter.second) {
+            inserter.first->second.push_back(globalinfo.polygon[*indexiter]);
+          }
+        }
+      }
+      // bonditer follows nicely
+      ASSERT( bonditer == bondlist.end(),
+          "FragmentationAction::performCall() - bonditer is out of step, it still points at bond "
+          +toString(*bonditer)+".");
+
+      // and insert
+      globalsaturationpositions.insert(
+          std::make_pair( _atom->getId(),
+              positions_per_neighbor
+          ));
+    }
+  }
 
   {