Changeset 4b2adf for src/Dynamics

Timestamp:

Jun 20, 2018, 8:20:43 AM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph_contraction-expansion, StoppableMakroAction

Children:

70caa3

Parents:

e3ce0e

git-author:

Frederik Heber <frederik.heber@…> (12/03/17 17:31:16)

git-committer:

Frederik Heber <frederik.heber@…> (06/20/18 08:20:43)

Message:

ForceAnnealiing now returns bool to indicate stop condition for MakroAction.

• Introducing FORCE_THRESHOLD below which we always stop.
• modified (simplified) bond side picking criterion.
• added some explanations to annealWithBondgraph().

File:

• src/Dynamics/ForceAnnealing.hpp (modified) (9 diffs)

src/Dynamics/ForceAnnealing.hpp

@@ -78 +78 @@
     maxSteps(_maxSteps),
     max_distance(_max_distance),
-    damping_factor(_damping_factor)
+    damping_factor(_damping_factor),
+    FORCE_THRESHOLD(1e-8)
   {}
 
@@ -94 +95 @@
    * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    * \param offset offset in matrix file to the first force component
+   * \return false - need to continue annealing, true - may stop because forces very small
    * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
    */
-  void operator()(
+  bool operator()(
       const int _TimeStep,
       const size_t _offset,
@@ -135 +137 @@
         << currentStep << " are " << maxComponents);
 
+    // check whether are smaller than threshold
+    bool AnnealingFinished = false;
+    double maxcomp = 0.;
+    for (size_t i=0;i<NDIM;++i)
+      maxcomp = std::max(maxcomp, fabs(maxComponents[i]));
+    if (maxcomp < FORCE_THRESHOLD) {
+      LOG(1, "STATUS: Force components are all less than " << FORCE_THRESHOLD
+          << ", stopping.");
+      currentStep = maxSteps;
+      AnnealingFinished = true;
+    }
+
     // are we in final step? Remember to reset static entities
     if (currentStep == maxSteps) {
@@ -140 +154 @@
       reset();
     }
+
+    return AnnealingFinished;
   }
 
@@ -431 +447 @@
       LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
 
+      /** for each atom, we imagine a plane at the position of the atom with
+       * its atomic gradient as the normal vector. We go through all its bonds
+       * and check on which side of the plane the bond is. This defines whether
+       * the bond is contracting (+) or expanding (-) with respect to this atom.
+       *
+       * A bond has two atoms, however. Hence, we do this for either atom and
+       * look at the combination: Is it in sum contracting or expanding given
+       * both projected_forces?
+       */
+
       /** go through all bonds and check projected_forces and side of plane
        * the idea is that if all bonds on one side are contracting ones or expanding,
@@ -440 +466 @@
        * if they are not mixed on either side, then we also only shift the
        * atom, namely away from expanding and towards contracting bonds.
+       *
+       * We may get this information right away by looking at the projected_forces.
+       * They give the atomic gradient of either atom projected onto the BondVector
+       * with an additional weight in [0,1].
        */
 
@@ -456 +486 @@
         std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
             projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside];
+        // note that projected_forces has sign such as to indicate whether
+        // atomic gradient wants bond to contract (-) or expand (+).
+        // This goes into sides: Minus side points away from gradient, plus side point
+        // towards gradient.
+        //
+        // the sum of both bond sides goes into types, depending on which is
+        // stronger if either wants a different thing
         const double &temp = forcelist[index];
         sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision
@@ -515 +552 @@
         LOG(4, "DEBUG: Chosen side is " << sideno << " with type " << typeside[sideno]);
         ASSERT( (typeside[sideno] == contracting) || (typeside[sideno] == expanding),
-            "annealWithBondGraph_BB() - chosen side is either expanding nor contracting.");
+            "annealWithBondGraph_BB() - chosen side is neither expanding nor contracting.");
         // one side is not mixed, all bonds on one side are of same type
         // hence, find out which bonds to exclude
         const BondList& ListOfBonds = walker.getListOfBonds();
 
-        // away from gradient (minus) and contracting
-        // or towards gradient (plus) and expanding
-        // gather all on same side and remove
-        const double sign =
-            (sides[first_per_side[sideno]] == types[first_per_side[sideno]])
-            ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]];
+        // sideno is away (0) or in direction (1) of gradient
+        // tpyes[first_per_side[sideno]] is either contracting (-1) or expanding (+1)
+        // : side (i), where (i) means which bonds we keep for the BFS, bonds
+        // on side (-i) are removed
+        // If all bonds on side away (0) want expansion (+1), move towards side with atom: side 1
+        // if all bonds side towards (1) want contraction (-1), move away side with atom : side -1
+
+        // unsure whether this or do nothing in the remaining cases:
+        // If all bonds on side toward (1) want expansion (+1), move away side with atom : side -1
+        //    (the reasoning is that the bond's other atom must have a stronger
+        //     gradient in the same direction and they push along atoms in
+        //     gradient direction: we don't want to interface with those.
+        //     Hence, move atoms along on away side
+        // if all bonds side away (0) want contraction (-1), move towards side with atom: side 1
+        //    (the reasoning is the same, don't interfere with update from
+        //     stronger gradient)
+        // hence, the decision is only based on sides once we have picked a side
+        // depending on all bonds associated with have same good type.
+        const double sign = sides[first_per_side[sideno]];
         LOG(4, "DEBUG: Removing edges from side with sign " << sign);
         BondList::const_iterator bonditer = ListOfBonds.begin();
@@ -649 +699 @@
   //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
   const double damping_factor;
+  //!> threshold for force components to stop annealing
+  const double FORCE_THRESHOLD;
 };