Index: src/Actions/MoleculeAction/ForceAnnealingAction.cpp
===================================================================
--- src/Actions/MoleculeAction/ForceAnnealingAction.cpp	(revision e3ce0edf75f86fbfc16237918f39d297d193fba2)
+++ src/Actions/MoleculeAction/ForceAnnealingAction.cpp	(revision 4b2adf09a4468da66e8b1d42eecb1e744c52878b)
@@ -35,4 +35,6 @@
 //#include "CodePatterns/MemDebug.hpp"
 
+#include "Actions/ActionExceptions.hpp"
+#include "Actions/MakroAction.hpp"
 #include "Actions/UndoRedoHelpers.hpp"
 #include "Atom/atom.hpp"
@@ -120,6 +122,20 @@
   // perform optimization step
   LOG(1, "Structural optimization.");
-  optimizer(CurrentStep, 1, params.UseBondGraph.get());
+  const bool StopStatus = optimizer(CurrentStep, 1, params.UseBondGraph.get());
   STATUS("Successfully optimized structure by one step.");
+
+  if (StopStatus && ActionQueue::getInstance().isMakroAction()) {
+    // send stop signal if we are taking part in MakroAction
+    MakroAction * const makroaction =
+        dynamic_cast<MakroAction *>(
+            const_cast<Action *>(
+                &ActionQueue::getInstance().getCurrentAction()));
+    if (makroaction != NULL) {
+      makroaction->setLoop(makroaction->getStep());
+    } else {
+      ELOG(2, "ActionQueue said we are inside process, but current Action is not a process?");
+      // do nothing
+    }
+  }
 
   std::vector< std::vector<AtomicInfo> > RedoInfo(2);
Index: src/Dynamics/ForceAnnealing.hpp
===================================================================
--- src/Dynamics/ForceAnnealing.hpp	(revision e3ce0edf75f86fbfc16237918f39d297d193fba2)
+++ src/Dynamics/ForceAnnealing.hpp	(revision 4b2adf09a4468da66e8b1d42eecb1e744c52878b)
@@ -78,5 +78,6 @@
     maxSteps(_maxSteps),
     max_distance(_max_distance),
-    damping_factor(_damping_factor)
+    damping_factor(_damping_factor),
+    FORCE_THRESHOLD(1e-8)
   {}
 
@@ -94,7 +95,8 @@
    * \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,4 +137,16 @@
         << 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,4 +154,6 @@
       reset();
     }
+
+    return AnnealingFinished;
   }
 
@@ -431,4 +447,14 @@
       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,4 +466,8 @@
        * 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,4 +486,11 @@
         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,15 +552,28 @@
         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,4 +699,6 @@
   //!> 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;
 };
 
