Changeset ce0ca4 for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

f5fa48

Parents:

c8d2e7

git-author:

Frederik Heber <heber@…> (07/18/14 16:24:53)

git-committer:

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

Message:

Added getAssociatedPoints().

• the general problem with getRemainingPoints() for saturating fragments with dangling bonds is that we violate the "Saturation Consistency Principle": Common saturation hydrogens for a specific fragments must remain at the exact same position for all containing fragments.
• only there do we ascertain that the eigenvalue is (due to the invariance of hydrogen to changes in its chemical neighborhood, valid to a good degree) actually removed by higher order fragments.
• this function is the first step to calculate a "global" set of saturation positions per atom.
• added also getIdentityAssociation() in case the association is trivial (i.e. in case of only bonds of degree 1).

File:

• src/Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (1 diff)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -1022 +1022 @@
 }
 
-
+SphericalPointDistribution::PolygonWithIndexTuples
+SphericalPointDistribution::getAssociatedPoints(
+    const WeightedPolygon_t &_polygon,
+    const int _N)
+{
+  SphericalPointDistribution::PolygonWithIndexTuples associatedpoints;
+
+  // initialze to given number of points
+  initSelf(_N);
+  LOG(2, "INFO: Matching old polygon " << _polygon
+      << " with new polygon " << points);
+
+  // check whether there are any points to associate
+  if (_polygon.empty()) {
+    associatedpoints.polygon.insert(
+        associatedpoints.polygon.end(),
+        points.begin(), points.end());
+    return associatedpoints;
+  }
+
+  if (_N > 0) {
+    // combine multiple points and create simple IndexList from IndexTupleList
+    MatchingControlStructure MCS = findBestMatching(_polygon);
+    IndexList_t bestmatching = joinPoints(points, MCS.newpoints, MCS.bestmatching);
+    LOG(2, "INFO: Best matching is " << bestmatching);
+
+    // gather the associated points (not the joined ones)
+    associatedpoints.polygon = MCS.newpoints;
+    // gather indices
+    associatedpoints.indices = MCS.bestmatching;
+
+    /// now we only need to rotate the associated points to match the given ones
+    /// with respect to the joined points in points
+
+    const size_t NumberIds = std::min(bestmatching.size(), (size_t)3);
+    // create old set
+    PolygonWithIndices oldSet;
+    oldSet.indices.resize(NumberIds, -1);
+    std::generate(oldSet.indices.begin(), oldSet.indices.end(), UniqueNumber);
+    for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
+        iter != _polygon.end(); ++iter)
+      oldSet.polygon.push_back(iter->first);
+
+    // _newpolygon has changed, so now convert to array with matched indices
+    PolygonWithIndices newSet;
+    SphericalPointDistribution::IndexList_t::const_iterator beginiter = bestmatching.begin();
+    SphericalPointDistribution::IndexList_t::const_iterator enditer = bestmatching.begin();
+    std::advance(enditer, NumberIds);
+    newSet.indices.resize(NumberIds, -1);
+    std::copy(beginiter, enditer, newSet.indices.begin());
+    std::copy(points.begin(),points.end(), std::back_inserter(newSet.polygon));
+
+    // determine rotation angles to align the two point distributions with
+    // respect to bestmatching:
+    // we use the center between the three first matching points
+    /// the first rotation brings these two centers to coincide
+    PolygonWithIndices rotatednewSet = newSet;
+    {
+      Rotation_t Rotation = findPlaneAligningRotation(oldSet, newSet);
+      LOG(5, "DEBUG: Rotating coordinate system by " << Rotation.second
+          << " around axis " << Rotation.first);
+      Line Axis(zeroVec, Rotation.first);
+
+      // apply rotation angle to bring newCenter to oldCenter in joined points
+      for (VectorArray_t::iterator iter = rotatednewSet.polygon.begin();
+          iter != rotatednewSet.polygon.end(); ++iter) {
+        Vector &current = *iter;
+        LOG(6, "DEBUG: Original joined point is " << current);
+        current =  Axis.rotateVector(current, Rotation.second);
+        LOG(6, "DEBUG: Rotated joined point is " << current);
+      }
+
+      // apply rotation angle to the whole set of associated points
+      for (VectorArray_t::iterator iter = associatedpoints.polygon.begin();
+          iter != associatedpoints.polygon.end(); ++iter) {
+        Vector &current = *iter;
+        LOG(6, "DEBUG: Original associated point is " << current);
+        current =  Axis.rotateVector(current, Rotation.second);
+        LOG(6, "DEBUG: Rotated associated point is " << current);
+      }
+
+#ifndef NDEBUG
+      // check: rotated "newCenter" should now equal oldCenter
+      // we don't check in case of two points as these lie on a great circle
+      // and the center cannot stably be recalculated. We may reactivate this
+      // when we calculate centers only once
+      if (oldSet.indices.size() > 2) {
+        Vector oldCenter;
+        Vector rotatednewCenter;
+        calculateOldAndNewCenters(
+            oldCenter, rotatednewCenter,
+            oldSet, rotatednewSet);
+        oldCenter.Normalize();
+        rotatednewCenter.Normalize();
+        // check whether centers are anti-parallel (factor -1)
+        // then we have the "non-unique poles" situation: points lie on great circle
+        // and both poles are valid solution
+        if (fabs(oldCenter.ScalarProduct(rotatednewCenter) + 1.)
+            < std::numeric_limits<double>::epsilon()*1e4)
+          rotatednewCenter *= -1.;
+        LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
+            << ", oldCenter is " << oldCenter);
+        const double difference = (rotatednewCenter - oldCenter).NormSquared();
+        ASSERT( difference < std::numeric_limits<double>::epsilon()*1e4,
+            "matchSphericalPointDistributions() - rotation does not work as expected by "
+            +toString(difference)+".");
+      }
+#endif
+    }
+    /// the second (orientation) rotation aligns the planes such that the
+    /// points themselves coincide
+    if (bestmatching.size() > 1) {
+      Rotation_t Rotation = findPointAligningRotation(oldSet, rotatednewSet);
+
+      // construct RotationAxis and two points on its plane, defining the angle
+      Rotation.first.Normalize();
+      const Line RotationAxis(zeroVec, Rotation.first);
+
+      LOG(5, "DEBUG: Rotating around self is " << Rotation.second
+          << " around axis " << RotationAxis);
+
+#ifndef NDEBUG
+      // check: first bestmatching in rotated_newpolygon and remainingnew
+      // should now equal
+      {
+        const IndexList_t::const_iterator iter = bestmatching.begin();
+
+        // check whether both old and newPosition are at same distance to oldCenter
+        Vector oldCenter = calculateCenter(oldSet);
+        const double distance = fabs(
+              (oldSet.polygon[0] - oldCenter).NormSquared()
+              - (rotatednewSet.polygon[*iter] - oldCenter).NormSquared()
+            );
+        LOG(4, "CHECK: Squared distance between oldPosition and newPosition "
+            << " with respect to oldCenter " << oldCenter << " is " << distance);
+//        ASSERT( distance < warn_amplitude,
+//            "matchSphericalPointDistributions() - old and newPosition's squared distance to oldCenter differs by "
+//            +toString(distance));
+
+        Vector rotatednew = RotationAxis.rotateVector(
+            rotatednewSet.polygon[*iter],
+            Rotation.second);
+        LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
+            << " while old was " << oldSet.polygon[0]);
+        const double difference = (rotatednew - oldSet.polygon[0]).NormSquared();
+        ASSERT( difference < distance+1e-8,
+            "matchSphericalPointDistributions() - orientation rotation ends up off by "
+            +toString(difference)+", more than "
+            +toString(distance+1e-8)+".");
+      }
+#endif
+
+      // align the set of associated points only here
+      for (VectorArray_t::iterator iter = associatedpoints.polygon.begin();
+          iter != associatedpoints.polygon.end(); ++iter) {
+        Vector &current = *iter;
+        LOG(6, "DEBUG: Original associated point is " << current);
+        current = RotationAxis.rotateVector(current, Rotation.second);
+        LOG(6, "DEBUG: Rotated associated point is " << current);
+      }
+    }
+  }
+
+  return associatedpoints;
+}
+
+SphericalPointDistribution::PolygonWithIndexTuples
+SphericalPointDistribution::getIdentityAssociation(
+    const WeightedPolygon_t &_polygon)
+{
+  unsigned int index = 0;
+  SphericalPointDistribution::PolygonWithIndexTuples returnpolygon;
+  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
+      iter != _polygon.end(); ++iter, ++index) {
+    returnpolygon.polygon.push_back( iter->first );
+    ASSERT( iter->second == 1,
+        "getIdentityAssociation() - bond with direction "
+        +toString(iter->second)
+        +" has degree higher than 1, getIdentityAssociation makes no sense.");
+    returnpolygon.indices.push_back( IndexList_t(1, index) );
+  }
+  return returnpolygon;
+}