Index: src/Dynamics/ForceAnnealing.hpp
===================================================================
--- src/Dynamics/ForceAnnealing.hpp	(revision 5f18f1e52e83873fccda63b356679d8150af3a2b)
+++ src/Dynamics/ForceAnnealing.hpp	(revision 4919f0f61d1ad1485a7fcee0cca520ea67762cf2)
@@ -92,17 +92,21 @@
    *
    *
-   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param _TimeStep current time step (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
    * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
    */
   void operator()(
-      const int _CurrentTimeStep,
+      const int _TimeStep,
       const size_t _offset,
       const bool _UseBondgraph)
   {
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT( CurrentTimeStep >= 0,
+        "ForceAnnealing::operator() - a new time step (upon which we work) must already have been copied.");
+
     // make sum of forces equal zero
     AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(
         _offset,
-        _CurrentTimeStep-1>=0 ? _CurrentTimeStep - 1 : 0);
+        CurrentTimeStep);
 
     // are we in initial step? Then set static entities
@@ -114,5 +118,5 @@
 
       // always use atomic annealing on first step
-      maxComponents = anneal(_CurrentTimeStep);
+      maxComponents = anneal(_TimeStep);
     } else {
       ++currentStep;
@@ -121,8 +125,8 @@
       // bond graph annealing is always followed by a normal annealing
       if (_UseBondgraph)
-        maxComponents = annealWithBondGraph(_CurrentTimeStep);
+        maxComponents = annealWithBondGraph_BarzilaiBorwein(_TimeStep);
       // cannot store RemnantGradient in Atom's Force as it ruins BB stepwidth calculation
       else
-        maxComponents = anneal(_CurrentTimeStep);
+        maxComponents = anneal_BarzilaiBorwein(_TimeStep);
     }
 
@@ -141,10 +145,16 @@
    * We assume that forces have just been calculated.
    *
-   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param _TimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    * \return to be filled with maximum force component over all atoms
    */
   Vector anneal(
-      const int CurrentTimeStep)
+      const int _TimeStep)
   {
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT( CurrentTimeStep >= 0,
+        "ForceAnnealing::anneal() - a new time step (upon which we work) must already have been copied.");
+
+    LOG(1, "STATUS: performing simple anneal with default stepwidth " << currentDeltat << " at step #" << currentStep);
+
     Vector maxComponents;
     bool deltat_decreased = false;
@@ -152,27 +162,12 @@
 		    iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
 			// atom's force vector gives steepest descent direction
-			const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
-			const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
-			const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
-			const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
-      LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
-      LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
-      LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
-      LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
+			const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
+			const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
+      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
+      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
 //			LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
 
       // we use Barzilai-Borwein update with position reversed to get descent
-      const Vector PositionDifference = currentPosition - oldPosition;
-			const Vector GradientDifference = (currentGradient - oldGradient);
-			double stepwidth = 0.;
-			if (GradientDifference.Norm() > MYEPSILON)
-        stepwidth = fabs(PositionDifference.ScalarProduct(GradientDifference))/
-            GradientDifference.NormSquared();
-			if (fabs(stepwidth) < 1e-10) {
-			  // dont' warn in first step, deltat usage normal
-			  if (currentStep != 1)
-          ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
-			  stepwidth = currentDeltat;
-			}
+			double stepwidth = currentDeltat;
       Vector PositionUpdate = stepwidth * currentGradient;
 			LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
@@ -186,5 +181,13 @@
 			// have different sign: Check whether this is the case and one step with
 			// deltat to interrupt this sequence
-			if ((currentStep > 1) && (!PositionDifference.IsZero()))
+			if (currentStep > 1) {
+		    const int OldTimeStep = _TimeStep-2;
+		    ASSERT( OldTimeStep >= 0,
+		        "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep)
+		        +", then there should be at least three time steps.");
+	      const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
+	      const Vector PositionDifference = currentPosition - oldPosition;
+        LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
+        LOG(4, "DEBUG: PositionDifference for atom #" << (*iter)->getId() << " is " << PositionDifference);
 				if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
 				    && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
@@ -198,9 +201,92 @@
 				      << ", using deltat: " << currentDeltat);
 				  PositionUpdate = currentDeltat * currentGradient;
+				}
 			}
 
 			// finally set new values
-			(*iter)->setPosition(currentPosition + PositionUpdate);
+			(*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
 		}
+
+    return maxComponents;
+  }
+
+  /** Performs Gradient optimization on the atoms using BarzilaiBorwein step width.
+   *
+   * \note this can only be called when there are at least two optimization
+   * time steps present, i.e. this must be preceeded by a simple anneal().
+   *
+   * We assume that forces have just been calculated.
+   *
+   * \param _TimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \return to be filled with maximum force component over all atoms
+   */
+  Vector anneal_BarzilaiBorwein(
+      const int _TimeStep)
+  {
+    const int OldTimeStep = _TimeStep-2;
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT( OldTimeStep >= 0,
+        "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
+    ASSERT(currentStep > 1,
+        "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
+
+    LOG(1, "STATUS: performing BarzilaiBorwein anneal at step #" << currentStep);
+
+    Vector maxComponents;
+    bool deltat_decreased = false;
+    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
+        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
+      // atom's force vector gives steepest descent direction
+      const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
+      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
+      const Vector oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
+      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
+      LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
+      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
+      LOG(4, "DEBUG: oldGradient for atom #" << (*iter)->getId() << " is " << oldGradient);
+      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
+//      LOG(4, "DEBUG: Force for atom #" << (*iter)->getId() << " is " << currentGradient);
+
+      // we use Barzilai-Borwein update with position reversed to get descent
+      const Vector PositionDifference = currentPosition - oldPosition;
+      const Vector GradientDifference = (currentGradient - oldGradient);
+      double stepwidth = 0.;
+      if (GradientDifference.Norm() > MYEPSILON)
+        stepwidth = fabs(PositionDifference.ScalarProduct(GradientDifference))/
+            GradientDifference.NormSquared();
+      if (fabs(stepwidth) < 1e-10) {
+        // dont' warn in first step, deltat usage normal
+        if (currentStep != 1)
+          ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
+        stepwidth = currentDeltat;
+      }
+      Vector PositionUpdate = stepwidth * currentGradient;
+      LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
+
+      // extract largest components for showing progress of annealing
+      for(size_t i=0;i<NDIM;++i)
+        if (currentGradient[i] > maxComponents[i])
+          maxComponents[i] = currentGradient[i];
+
+      // steps may go back and forth again (updates are of same magnitude but
+      // have different sign: Check whether this is the case and one step with
+      // deltat to interrupt this sequence
+      if (!PositionDifference.IsZero())
+        if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
+            && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
+          // for convergence we want a null sequence here, too
+          if (!deltat_decreased) {
+            deltat_decreased = true;
+            currentDeltat = .5*currentDeltat;
+          }
+          LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate
+              << " > " << PositionDifference
+              << ", using deltat: " << currentDeltat);
+          PositionUpdate = currentDeltat * currentGradient;
+      }
+
+      // finally set new values
+      (*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
+    }
 
     return maxComponents;
@@ -244,5 +330,8 @@
   }
 
-  /** Performs Gradient optimization on the bonds.
+  /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
+   *
+   * \note this can only be called when there are at least two optimization
+   * time steps present, i.e. this must be preceeded by a simple anneal().
    *
    * We assume that forces have just been calculated. These forces are projected
@@ -251,10 +340,19 @@
    *
    *
-   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param _TimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    * \param maxComponents to be filled with maximum force component over all atoms
    */
-  Vector annealWithBondGraph(
-      const int CurrentTimeStep)
+  Vector annealWithBondGraph_BarzilaiBorwein(
+      const int _TimeStep)
   {
+    const int OldTimeStep = _TimeStep-2;
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT(currentStep > 1,
+        "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
+    ASSERT(OldTimeStep >= 0,
+        "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth, and the new one to update on already present.");
+
+    LOG(1, "STATUS: performing BarzilaiBorwein anneal on bonds at step #" << currentStep);
+
     Vector maxComponents;
 
@@ -297,5 +395,5 @@
         AtomicForceManipulator<T>::atoms.begin(),
         AtomicForceManipulator<T>::atoms.end(),
-        CurrentTimeStep);
+        _TimeStep);
     const BondVectors::container_t &sorted_bonds = bv.getSorted();
 
@@ -316,13 +414,13 @@
     RemnantGradient_per_atom_t RemnantGradient_per_atom;
     for (size_t timestep = 0; timestep <= 1; ++timestep) {
-      const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
-      LOG(2, "DEBUG: CurrentTimeStep is " << CurrentTimeStep
+      const size_t ReferenceTimeStep = _TimeStep-timestep-1;
+      LOG(2, "DEBUG: given time step is " << _TimeStep
           << ", timestep is " << timestep
-          << ", and CurrentStep is " << CurrentStep);
+          << ", and ReferenceTimeStep is " << ReferenceTimeStep);
 
       for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
           iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
         const atom &walker = *(*iter);
-        const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
+        const Vector &walkerGradient = walker.getAtomicForceAtStep(ReferenceTimeStep);
         LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
             << walker << " is " << walkerGradient);
@@ -333,5 +431,5 @@
           // gather subset of BondVectors for the current atom
           const std::vector<Vector> BondVectors =
-              bv.getAtomsBondVectorsAtStep(walker, CurrentStep);
+              bv.getAtomsBondVectorsAtStep(walker, ReferenceTimeStep);
 
           // go through all its bonds and calculate what magnitude is represented
@@ -340,5 +438,5 @@
               weights_per_atom[timestep].insert(
                   std::make_pair(walker.getId(),
-                  bv.getWeightsForAtomAtStep(walker, BondVectors, CurrentStep)) );
+                  bv.getWeightsForAtomAtStep(walker, BondVectors, ReferenceTimeStep)) );
           ASSERT( inserter.second,
               "ForceAnnealing::operator() - weight map for atom "+toString(walker)
@@ -428,6 +526,6 @@
         LOG(4, "DEBUG: current projected gradient for "
             << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
-        const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
-        const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
+        const Vector &oldPosition = bondatom[side]->getPositionAtStep(OldTimeStep);
+        const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep);
         const Vector PositionDifference = currentPosition - oldPosition;
         LOG(4, "DEBUG: old position is " << oldPosition);
@@ -478,5 +576,5 @@
       atom &walker = *(*iter);
       // extract largest components for showing progress of annealing
-      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
+      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
       for(size_t i=0;i<NDIM;++i)
         if (currentGradient[i] > maxComponents[i])
@@ -505,7 +603,6 @@
       LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
           << ", namely " << *walker);
-      walker->setPosition(
-          walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
-          + update - CommonTranslation);
+      walker->setPositionAtStep(_TimeStep,
+          walker->getPositionAtStep(CurrentTimeStep) + update - CommonTranslation);
 //      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
     }
