Changeset 08df4a for src/Dynamics/ForceAnnealing.hpp

Timestamp:

Aug 14, 2017, 8:12:49 AM (8 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_with_BondGraph_continued

Children:

a6574a

Parents:

0cc08c

git-author:

Frederik Heber <frederik.heber@…> (08/10/17 15:35:45)

git-committer:

Frederik Heber <frederik.heber@…> (08/14/17 08:12:49)

Message:

Rewrote annealWithBondGraph_BarzilaiBorwein() to simply distinguish expansion and contraction in bonds.

• we shift the neighboring set away in case of expansion and towards in case of contraction.

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -185 +185 @@
                         "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep)
                         +", then there should be at least three time steps.");
-              const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
+              const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
               const Vector PositionDifference = currentPosition - oldPosition;
         LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
@@ -237 +237 @@
         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);
+      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);
@@ -291 +291 @@
   }
 
-  // knowing the number of bonds in total, we can setup the storage for the
-  // projected forces
-  enum whichatom_t {
-    leftside=0,
-    rightside=1,
-    MAX_sides
-  };
-
-  /** Helper function to put bond force into a container.
-   *
-   * \param _walker atom
-   * \param _current_bond current bond of \a _walker
-   * \param _timestep time step
-   * \param _force calculated bond force
-   * \param _bv bondvectors for obtaining the correct index
-   * \param _projected_forces container
-   */
-  static void ForceStore(
-      const atom &_walker,
-      const bond::ptr &_current_bond,
-      const size_t &_timestep,
-      const double _force,
-      const BondVectors &_bv,
-      std::vector< // time step
-            std::vector< // which bond side
-              std::vector<double> > // over all bonds
-                > &_projected_forces)
-  {
-    std::vector<double> &forcelist = (&_walker == _current_bond->leftatom) ?
-        _projected_forces[_timestep][leftside] : _projected_forces[_timestep][rightside];
-    const size_t index = _bv.getIndexForBond(_current_bond);
-    ASSERT( index != (size_t)-1,
-        "ForceAnnealing() - could not find bond "+toString(*_current_bond)
-        +" in bondvectors");
-    forcelist[index] = _force;
-  }
-
   /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
    *
@@ -364 +327 @@
     BreadthFirstSearchGatherer NodeGatherer(BGcreator);
 
-    /// We assume that a force is local, i.e. a bond is too short yet and hence
-    /// the atom needs to be moved. However, all the adjacent (bound) atoms might
-    /// already be at the perfect distance. If we just move the atom alone, we ruin
-    /// all the other bonds. Hence, it would be sensible to move every atom found
-    /// through the bond graph in the direction of the force as well by the same
-    /// PositionUpdate. This is almost what we are going to do.
-
-    /// One issue is: If we need to shorten bond, then we use the PositionUpdate
-    /// also on the the other bond partner already. This is because it is in the
-    /// direction of the bond. Therefore, the update is actually performed twice on
-    /// each bond partner, i.e. the step size is twice as large as it should be.
-    /// This problem only occurs when bonds need to be shortened, not when they
-    /// need to be made longer (then the force vector is facing the other
-    /// direction than the bond vector).
-    /// As a remedy we need to average the force on either end of the bond and
-    /// check whether each gradient points inwards out or outwards with respect
-    /// to the bond and then shift accordingly.
-
-    /// One more issue is that the projection onto the bond directions does not
-    /// recover the gradient but may be larger as the bond directions are a
-    /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would
-    /// suffice to span the plane). To this end, we need to account for the
-    /// overestimation and obtain a weighting for each bond.
+    /** We assume that a force is local, i.e. a bond is too short yet and hence
+     * the atom needs to be moved. However, all the adjacent (bound) atoms might
+     * already be at the perfect distance. If we just move the atom alone, we ruin
+     * all the other bonds. Hence, it would be sensible to move every atom found
+     * through the bond graph in the direction of the force as well by the same
+     * PositionUpdate. This is almost what we are going to do, see below.
+     *
+     * This is to make the force a little more global in the sense of a multigrid
+     * solver that uses various coarser grids to transport errors more effectively
+     * over finely resolved grids.
+     *
+     */
+
+    /** The idea is that we project the gradients onto the bond vectors and determine
+     * from the sum of projected gradients from either side whether the bond is
+     * to contract or to expand. As the gradient acting as the normal vector of
+     * a plane supported at the position of the atom separates all bonds into two
+     * sets, we check whether all on one side are contracting and all on the other
+     * side are expanding. In this case we may move not only the atom itself but
+     * may propagate its update along a limited-horizon BFS to neighboring atoms.
+     *
+     */
 
     // initialize helper class for bond vectors using bonds from range of atoms
@@ -394 +356 @@
         AtomicForceManipulator<T>::atoms.end(),
         _TimeStep);
-    const BondVectors::container_t &sorted_bonds = bv.getSorted();
-
-    std::vector< // time step
-      std::vector< // which bond side
-        std::vector<double> > // over all bonds
-          > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps)
-    for (size_t i=0;i<2;++i) {
-      projected_forces[i].resize(MAX_sides);
-      for (size_t j=0;j<MAX_sides;++j)
-        projected_forces[i][j].resize(sorted_bonds.size(), 0.);
+
+    std::vector< // which bond side
+      std::vector<double> > // over all bonds
+        projected_forces; // one for leftatoms, one for rightatoms
+    projected_forces.resize(BondVectors::MAX_sides);
+    for (size_t j=0;j<BondVectors::MAX_sides;++j)
+      projected_forces[j].resize(bv.size(), 0.);
+
+    // for each atom we need to project the gradient
+    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(CurrentTimeStep);
+      const double GradientNorm = walkerGradient.Norm();
+      LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
+          << walker << " is " << walkerGradient << " with magnitude of "
+          << GradientNorm);
+
+      if (GradientNorm > MYEPSILON) {
+        bv.getProjectedGradientsForAtomAtStep(
+            walker, walkerGradient, CurrentTimeStep, projected_forces
+        );
+      } else {
+        LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
+            << MYEPSILON << " for atom " << walker);
+        // note that projected_forces is initialized to full length and filled
+        // with zeros. Hence, nothing to do here
+      }
     }
 
-    // for each atom we need to gather weights and then project the gradient
-    typedef std::map<atomId_t, BondVectors::weights_t > weights_per_atom_t;
-    std::vector<weights_per_atom_t> weights_per_atom(2);
-    typedef std::map<atomId_t, Vector> RemnantGradient_per_atom_t;
-    RemnantGradient_per_atom_t RemnantGradient_per_atom;
-    for (size_t timestep = 0; timestep <= 1; ++timestep) {
-      const size_t ReferenceTimeStep = _TimeStep-timestep-1;
-      LOG(2, "DEBUG: given time step is " << _TimeStep
-          << ", timestep is " << timestep
-          << ", 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(ReferenceTimeStep);
-        LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
-            << walker << " is " << walkerGradient << " with magnitude of "
-            << walkerGradient.Norm());
-
+    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
+    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
+        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
+      atom &walker = *(*iter);
+
+      /// calculate step width
+      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);
+
+      /** 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,
+       * respectively, then we may shift not only the atom with respect to its
+       * gradient but also its neighbors (towards contraction or towards
+       * expansion depending on direction of gradient).
+       * if they are mixed on both sides of the plane, then we simply shift
+       * only the atom itself.
+       * if they are not mixed on either side, then we also only shift the
+       * atom, namely away from expanding and towards contracting bonds.
+       */
+
+      // sign encodes side of plane and also encodes contracting(-) or expanding(+)
+      typedef std::vector<int> sides_t;
+      typedef std::vector<int> types_t;
+      sides_t sides;
+      types_t types;
+      const BondList& ListOfBonds = walker.getListOfBonds();
+      for(BondList::const_iterator bonditer = ListOfBonds.begin();
+          bonditer != ListOfBonds.end(); ++bonditer) {
+        const bond::ptr &current_bond = *bonditer;
+
+        // BondVector goes from bond::rightatom to bond::leftatom
+        const size_t index = bv.getIndexForBond(current_bond);
+        std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
+            projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside];
+        const double &temp = forcelist[index];
+        sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision
+        const double sum =
+            projected_forces[BondVectors::leftside][index]+projected_forces[BondVectors::rightside][index];
+        types.push_back( sum/fabs(sum) );
+        LOG(4, "DEBUG: Bond " << *current_bond << " is on side " << sides.back()
+            << " and has type " << types.back());
+      }
+//      /// check whether both conditions are compatible:
+//      // i.e. either we have ++/-- for all entries in sides and types
+//      // or we have +-/-+ for all entries
+//      // hence, multiplying and taking the sum and its absolute value
+//      // should be equal to the maximum number of entries
+//      sides_t results;
+//      std::transform(
+//          sides.begin(), sides.end(),
+//          types.begin(),
+//          std::back_inserter(results),
+//          std::multiplies<int>);
+//      int result = abs(std::accumulate(results.begin(), results.end(), 0, std::plus<int>));
+
+      std::vector<size_t> first_per_side(2, (size_t)-1); //!< mark down one representative from either side
+      std::vector< std::vector<int> > types_per_side(2); //!< gather all types on each side
+      types_t::const_iterator typesiter = types.begin();
+      for (sides_t::const_iterator sidesiter = sides.begin();
+          sidesiter != sides.end(); ++sidesiter, ++typesiter) {
+        const size_t index = (*sidesiter+1)/2;
+        types_per_side[index].push_back(*typesiter);
+        if (first_per_side[index] == (size_t)-1)
+          first_per_side[index] = std::distance(const_cast<const sides_t &>(sides).begin(), sidesiter);
+      }
+      LOG(4, "DEBUG: First on side minus is " << first_per_side[0] << ", and first on side plus is "
+          << first_per_side[1]);
+      //!> enumerate types per side with a little witching with the numbers to allow easy setting from types
+      enum whichtypes_t {
+        contracting=0,
+        unset=1,
+        expanding=2,
+        mixed,
+        MAX_types
+      };
+      std::vector<int> typeside(2, unset);
+      for(size_t i=0;i<2;++i) {
+        for (std::vector<int>::const_iterator tpsiter = types_per_side[i].begin();
+            tpsiter != types_per_side[i].end(); ++tpsiter) {
+          if (typeside[i] == unset) {
+            typeside[i] = *tpsiter+1; //contracting(0) or expanding(2)
+          } else {
+            if (typeside[i] != (*tpsiter+1)) // no longer he same type
+              typeside[i] = mixed;
+          }
+        }
+      }
+      LOG(4, "DEBUG: Minus side is " << typeside[0] << " and plus side is " << typeside[1]);
+
+      // check whether positive side is contraction, then pick negative side
+      int sideno = 1;
+      if (typeside[1] == contracting) {
+        sideno = 0;
+      }
+
+      typedef std::vector< std::pair<atomId_t, atomId_t> > RemovedEdges_t;
+      RemovedEdges_t RemovedEdges;
+      if (sideno >= 0) {
+        sideno = ((sideno+1) % 2)*2-1; // invert and make it -1/+1
         const BondList& ListOfBonds = walker.getListOfBonds();
-        if (walkerGradient.Norm() > MYEPSILON) {
-
-          // gather subset of BondVectors for the current atom
-          const std::vector<Vector> BondVectors =
-              bv.getAtomsBondVectorsAtStep(walker, ReferenceTimeStep);
-
-          // go through all its bonds and calculate what magnitude is represented
-          // by the others i.e. sum of scalar products against other bonds
-          const std::pair<weights_per_atom_t::iterator, bool> inserter =
-              weights_per_atom[timestep].insert(
-                  std::make_pair(walker.getId(),
-                  bv.getWeightsForAtomAtStep(walker, BondVectors, ReferenceTimeStep)) );
-          ASSERT( inserter.second,
-              "ForceAnnealing::operator() - weight map for atom "+toString(walker)
-              +" and time step "+toString(timestep)+" already filled?");
-          BondVectors::weights_t &weights = inserter.first->second;
-          ASSERT( weights.size() == ListOfBonds.size(),
-              "ForceAnnealing::operator() - number of weights "
-              +toString(weights.size())+" does not match number of bonds "
-              +toString(ListOfBonds.size())+", error in calculation?");
-
-          // projected gradient over all bonds and place in one of projected_forces
-          // using the obtained weights
-          BondVectors::forcestore_t forcestoring =
-              boost::bind(&ForceAnnealing::ForceStore, _1, _2, _3, _4,
-                  boost::cref(bv), boost::ref(projected_forces));
-          const Vector RemnantGradient = bv.getRemnantGradientForAtomAtStep(
-              walker, walkerGradient, BondVectors, weights, timestep, forcestoring
-          );
-          RemnantGradient_per_atom.insert( std::make_pair(walker.getId(), RemnantGradient) );
-        } else {
-          LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
-              << MYEPSILON << " for atom " << walker);
-          // note that projected_forces is initialized to full length and filled
-          // with zeros. Hence, nothing to do here
+
+        BondList::const_iterator bonditer = ListOfBonds.begin();
+        for (sides_t::const_iterator sidesiter = sides.begin();
+            sidesiter != sides.end(); ++sidesiter, ++bonditer) {
+          if (*sidesiter == sideno) {
+            // remove the edge
+            const bond::ptr &current_bond = *bonditer;
+            RemovedEdges.push_back( std::make_pair(
+                current_bond->leftatom->getId(),
+                current_bond->rightatom->getId())
+            );
+            LOG(5, "DEBUG: Removing edge " << RemovedEdges.back());
+#ifndef NDEBUG
+            const bool status =
+#endif
+                BGcreator.removeEdge(RemovedEdges.back());
+            ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
+          }
         }
-      }
-    }
-
-    // step through each bond and shift the atoms
-    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
-
-    LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
-    const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep);
-
-    for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
-        bondsiter != sorted_bonds.end(); ++bondsiter) {
-      const bond::ptr &current_bond = *bondsiter;
-      const size_t index = bv.getIndexForBond(current_bond);
-      const atom* bondatom[MAX_sides] = {
-          current_bond->leftatom,
-          current_bond->rightatom
-      };
-
-      // remove the edge
-#ifndef NDEBUG
-      const bool status =
-#endif
-          BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
-      ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
-
-      // gather nodes for either atom
-      BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
-      BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
-      for (size_t side=leftside;side<MAX_sides;++side) {
-        bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
-        distance_map[side] = NodeGatherer.getDistances();
-        std::sort(bondside_set[side].begin(), bondside_set[side].end());
-      }
-
-      // re-add edge
-      BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
-
-      // do for both leftatom and rightatom of bond
-      for (size_t side = leftside; side < MAX_sides; ++side) {
-        const double &bondforce = projected_forces[0][side][index];
-        const double &oldbondforce = projected_forces[1][side][index];
-        const double bondforcedifference = fabs(bondforce - oldbondforce);
-        LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
-            << " side of bond is " << bondforce);
-        LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
-            << " side of bond is " << oldbondforce);
-        // if difference or bondforce itself is zero, do nothing
-        if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
-          continue;
-
-        // get BondVector to bond
-        const BondVectors::mapped_t::const_iterator bviter =
-            bondvectors.find(current_bond);
-        ASSERT( bviter != bondvectors.end(),
-            "ForceAnnealing() - cannot find current_bond ?");
-        ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
-            "ForceAnnealing() - norm of BondVector is not one");
-        const Vector &BondVector = bviter->second;
-
-        // calculate gradient and position differences for stepwidth
-        const Vector currentGradient = bondforce * BondVector;
-        LOG(4, "DEBUG: current projected gradient for "
-            << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
-        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);
-        LOG(4, "DEBUG: current position is " << currentPosition);
-        LOG(4, "DEBUG: difference in position is " << PositionDifference);
-        LOG(4, "DEBUG: bondvector is " << BondVector);
-        const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector);
-        LOG(4, "DEBUG: difference in position projected onto bondvector is "
-            << projected_PositionDifference);
-        LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
-
-        // calculate step width
-        double stepwidth =
-            fabs(projected_PositionDifference)/bondforcedifference;
-        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);
-
-        // add PositionUpdate for all nodes in the bondside_set
-        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
-            setiter != bondside_set[side].end(); ++setiter) {
+        // perform limited-horizon BFS
+        BoostGraphHelpers::Nodeset_t bondside_set;
+        BreadthFirstSearchGatherer::distance_map_t distance_map;
+        bondside_set = NodeGatherer(walker.getId(), max_distance);
+        distance_map = NodeGatherer.getDistances();
+        std::sort(bondside_set.begin(), bondside_set.end());
+
+        // re-add edges
+        for (RemovedEdges_t::const_iterator iter = RemovedEdges.begin();
+            iter != RemovedEdges.end(); ++iter)
+          BGcreator.addEdge(*iter);
+        RemovedEdges.clear();
+
+        // update position with dampening factor on the discovered bonds
+        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
+            setiter != bondside_set.end(); ++setiter) {
           const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
-            = distance_map[side].find(*setiter);
-          ASSERT( diter != distance_map[side].end(),
+            = distance_map.find(*setiter);
+          ASSERT( diter != distance_map.end(),
               "ForceAnnealing() - could not find distance to an atom.");
           const double factor = pow(damping_factor, diter->second+1);
@@ -568 +560 @@
           }
         }
+      } else {
+        // simple atomic annealing, i.e. damping factor of 1
+        LOG(3, "DEBUG: Update for atom #" << walker.getId() << " will be " << PositionUpdate);
+        GatheredUpdates.insert(
+            std::make_pair(
+                walker.getId(),
+                PositionUpdate) );
       }
     }
@@ -580 +579 @@
     }
 
-    // remove center of weight translation from gathered updates
-    Vector CommonTranslation;
-    for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
-        iter != GatheredUpdates.end(); ++iter) {
-      const Vector &update = iter->second;
-      CommonTranslation += update;
-    }
-    CommonTranslation *= 1./(double)GatheredUpdates.size();
-    LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation
-        << " from all updates.");
+//    // remove center of weight translation from gathered updates
+//    Vector CommonTranslation;
+//    for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
+//        iter != GatheredUpdates.end(); ++iter) {
+//      const Vector &update = iter->second;
+//      CommonTranslation += update;
+//    }
+//    CommonTranslation *= 1./(double)GatheredUpdates.size();
+//    LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation
+//        << " from all updates.");
 
     // apply the gathered updates and set remnant gradients for atomic annealing
@@ -602 +601 @@
           << ", namely " << *walker);
       walker->setPositionAtStep(_TimeStep,
-          walker->getPositionAtStep(CurrentTimeStep) + update - CommonTranslation);
-//      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
+          walker->getPositionAtStep(CurrentTimeStep) + update); // - CommonTranslation);
     }