Changeset 8754b1 for src/Dynamics/ForceAnnealing.hpp

Timestamp:

Jul 5, 2017, 7:46:00 PM (8 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_oldresults

Parents:

bd0159

git-author:

Frederik Heber <frederik.heber@…> (05/19/17 14:29:40)

git-committer:

Frederik Heber <frederik.heber@…> (07/05/17 19:46:00)

Message:

ForceAnnealing now shifts atoms in the bond graph to anneal forces, too.

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -20 +20 @@
 #include "CodePatterns/Log.hpp"
 #include "CodePatterns/Verbose.hpp"
+#include "Descriptors/AtomIdDescriptor.hpp"
 #include "Dynamics/AtomicForceManipulator.hpp"
 #include "Fragmentation/ForceMatrix.hpp"
+#include "Graph/BoostGraphCreator.hpp"
+#include "Graph/BoostGraphHelpers.hpp"
+#include "Graph/BreadthFirstSearchGatherer.hpp"
 #include "Helpers/helpers.hpp"
 #include "Helpers/defs.hpp"
@@ -29 +33 @@
 #include "World.hpp"
 
-/** This class is the essential build block for performing structual optimization.
+/** This class is the essential build block for performing structural optimization.
  *
  * Sadly, we have to use some static instances as so far values cannot be passed
  * between actions. Hence, we need to store the current step and the adaptive-
- * step width (we cannot perform a linesearch, as we have no control over the
+ * step width (we cannot perform a line search, as we have no control over the
  * calculation of the forces).
+ *
+ * However, we do use the bond graph, i.e. if a single atom needs to be shifted
+ * to the left, then the whole molecule left of it is shifted, too. This is
+ * controlled by the \a max_distance parameter.
  */
 template <class T>
@@ -48 +56 @@
    * \param _IsAngstroem whether length units are in angstroem or bohr radii
    * \param _maxSteps number of optimization steps to perform
+   * \param _max_distance up to this bond order is bond graph taken into account.
    */
   ForceAnnealing(
       AtomSetMixin<T> &_atoms,
       bool _IsAngstroem,
-      const size_t _maxSteps) :
+      const size_t _maxSteps,
+      const int _max_distance) :
     AtomicForceManipulator<T>(_atoms, 1., _IsAngstroem),
-    maxSteps(_maxSteps)
+    maxSteps(_maxSteps),
+    max_distance(_max_distance),
+    damping_factor(0.5)
   {}
   /** Destructor of class ForceAnnealing.
@@ -86 +98 @@
     }
 
+    // get nodes on either side of selected bond via BFS discovery
+//    std::vector<atomId_t> atomids;
+//    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
+//        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
+//      atomids.push_back((*iter)->getId());
+//    }
+//    ASSERT( atomids.size() == AtomicForceManipulator<T>::atoms.size(),
+//        "ForceAnnealing() - could not gather all atomic ids?");
+    BoostGraphCreator BGcreator;
+    BGcreator.createFromRange(
+        AtomicForceManipulator<T>::atoms.begin(),
+        AtomicForceManipulator<T>::atoms.end(),
+        AtomicForceManipulator<T>::atoms.size(),
+        BreadthFirstSearchGatherer::AlwaysTruePredicate);
+    BreadthFirstSearchGatherer NodeGatherer(BGcreator);
+
+    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
     Vector maxComponents(zeroVec);
     for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
@@ -110 +139 @@
       LOG(3, "DEBUG: Update would be " << PositionUpdate);
 
+//      // add update to central atom
+//      const atomId_t atomid = (*iter)->getId();
+//      if (GatheredUpdates.count(atomid)) {
+//        GatheredUpdates[atomid] += PositionUpdate;
+//      } else
+//        GatheredUpdates.insert( std::make_pair(atomid, PositionUpdate) );
+
+      // 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 just what we are going to do.
+
+      /// get all nodes from bonds in the direction of the current force
+
+      // remove edges facing in the wrong direction
+      std::vector<bond::ptr> removed_bonds;
+      const BondList& ListOfBonds = (*iter)->getListOfBonds();
+      for(BondList::const_iterator bonditer = ListOfBonds.begin();
+          bonditer != ListOfBonds.end(); ++bonditer) {
+        const bond &current_bond = *(*bonditer);
+        LOG(2, "DEBUG: Looking at bond " << current_bond);
+        Vector BondVector = (*iter)->getPosition();
+        BondVector -= ((*iter)->getId() == current_bond.rightatom->getId())
+            ? current_bond.rightatom->getPosition() : current_bond.leftatom->getPosition();
+        BondVector.Normalize();
+        if (BondVector.ScalarProduct(currentGradient) < 0) {
+          removed_bonds.push_back(*bonditer);
+#ifndef NDEBUG
+          const bool status =
+#endif
+              BGcreator.removeEdge(current_bond.leftatom->getId(), current_bond.rightatom->getId());
+          ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
+        }
+      }
+      BoostGraphHelpers::Nodeset_t bondside_set = NodeGatherer((*iter)->getId(), max_distance);
+      const BreadthFirstSearchGatherer::distance_map_t &distance_map = NodeGatherer.getDistances();
+      std::sort(bondside_set.begin(), bondside_set.end());
+      // re-add those edges
+      for (std::vector<bond::ptr>::const_iterator bonditer = removed_bonds.begin();
+          bonditer != removed_bonds.end(); ++bonditer)
+        BGcreator.addEdge((*bonditer)->leftatom->getId(), (*bonditer)->rightatom->getId());
+
+      // apply PositionUpdate to all nodes in the bondside_set
+      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.find(*setiter);
+        ASSERT( diter != distance_map.end(),
+            "ForceAnnealing() - could not find distance to an atom.");
+        const double factor = pow(damping_factor, diter->second);
+        LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
+            << factor << "*" << PositionUpdate);
+        if (GatheredUpdates.count((*setiter))) {
+          GatheredUpdates[(*setiter)] += factor*PositionUpdate;
+        } else {
+          GatheredUpdates.insert(
+              std::make_pair(
+                  (*setiter),
+                  factor*PositionUpdate) );
+        }
+      }
+
       // extract largest components for showing progress of annealing
       for(size_t i=0;i<NDIM;++i)
@@ -115 +208 @@
           maxComponents[i] = currentGradient[i];
 
-      // are we in initial step? Then don't check against velocity
-      if ((currentStep > 1) && (!(*iter)->getAtomicVelocity().IsZero()))
-        // update with currentDelta tells us how the current gradient relates to
-        // the last one: If it has become larger, reduce currentDelta
-        if ((PositionUpdate.ScalarProduct((*iter)->getAtomicVelocity()) < 0)
-            && (currentDeltat > MinimumDeltat)) {
-          currentDeltat = .5*currentDeltat;
-          LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate.NormSquared()
-              << " > " << (*iter)->getAtomicVelocity().NormSquared()
-              << ", decreasing deltat: " << currentDeltat);
-          PositionUpdate = currentDeltat * currentGradient;
-        }
-
-      // finally set new values
-      (*iter)->setPosition(currentPosition + PositionUpdate);
-      (*iter)->setAtomicVelocity(PositionUpdate);
-      //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
-//        (*iter)->VelocityVerletUpdateU((*iter)->getId(), NextStep-1, Deltat, IsAngstroem);
-
       // reset force vector for next step except on final one
       if (currentStep != maxSteps)
         (*iter)->setAtomicForce(zeroVec);
+    }
+    // apply the gathered updates
+    for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
+        iter != GatheredUpdates.end(); ++iter) {
+      const atomId_t &atomid = iter->first;
+      const Vector &update = iter->second;
+      atom* const walker = World::getInstance().getAtom(AtomById(atomid));
+      ASSERT( walker != NULL,
+          "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
+      walker->setPosition( walker->getPosition() + update );
+      walker->setAtomicVelocity(update);
+      LOG(3, "DEBUG: Applying update " << update << " to atom " << *walker);
     }
 
@@ -171 +257 @@
   //!> minimum deltat for internal while loop (adaptive step width)
   static double MinimumDeltat;
+  //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
+  const int max_distance;
+  //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
+  const double damping_factor;
 };