Changeset 8754b1

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.

Files:

• src/Actions/MoleculeAction/ForceAnnealingAction.cpp (modified) (1 diff)
• src/Actions/MoleculeAction/ForceAnnealingAction.def (modified) (1 diff)
• src/Dynamics/ForceAnnealing.hpp (modified) (7 diffs)
• tests/Python/AllActions/options.dat (modified) (1 diff)

src/Actions/MoleculeAction/ForceAnnealingAction.cpp

@@ -93 +93 @@
       set,
       true,
-      params.steps.get());
+      params.steps.get(),
+      params.MaxDistance.get());
   // parse forces into next step
   if (!params.forcesfile.get().string().empty()) {

src/Actions/MoleculeAction/ForceAnnealingAction.def

@@ -16 +16 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
-#define paramtypes (boost::filesystem::path)(unsigned int)(bool)
-#define paramtokens ("forces-file")("steps")("output-every-step")
-#define paramdescriptions ("file containing")("fixed number of optimization steps to be performed")("whether WorldTime should be increased and output written after every step, useful if optimization might hang")
-#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(PARAM_DEFAULT("0"))
-#define paramreferences (forcesfile)(steps)(DoOutput)
+#define paramtypes (boost::filesystem::path)(unsigned int)(bool)(int)
+#define paramtokens ("forces-file")("steps")("output-every-step")("max-distance")
+#define paramdescriptions ("file containing")("fixed number of optimization steps to be performed")("whether WorldTime should be increased and output written after every step, useful if optimization might hang")("maximum distance to which bond graph is taken into account")
+#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(PARAM_DEFAULT("0"))(PARAM_DEFAULT(0))
+#define paramreferences (forcesfile)(steps)(DoOutput)(MaxDistance)
 #define paramvalids \
 (DummyValidator< boost::filesystem::path >()) \
 (NotZeroValidator< unsigned int >()) \
-(DummyValidator<bool>())
+(DummyValidator<bool>()) \
+(DummyValidator< int >())
 
 #define statetypes (std::vector<AtomicInfo>)(std::vector<AtomicInfo>)

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;
 };
 

tests/Python/AllActions/options.dat

@@ -104 +104 @@
 load-session    "test.py"
 MaxDistance     "-1"
+max-distance    "0"
 max-meshwidth   "0.3"
 mesh-offset     "0.5,0.5,0.5"