Changeset 6d360a for src

Timestamp:

Jul 20, 2017, 9:38:38 AM (8 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_with_BondGraph_continued

Children:

9d3846

Parents:

4cd46b

git-author:

Frederik Heber <frederik.heber@…> (06/20/17 20:26:12)

git-committer:

Frederik Heber <frederik.heber@…> (07/20/17 09:38:38)

Message:

ForceAnnealing can now be used either atom- or bond-centered.

• new keyword "use-bondgraph" to use either case.
• atom is old annealing method, bond graph is new annealing method where the bond neighborhood is shifted as well to anneal to force projected onto the BondVector.
• In the first step we always use atom-centered annealing where we use the given deltat step width. Afterwards Barzilai-Borwein may then proceed.

Location:

src

Files:

• Actions/MoleculeAction/ForceAnnealingAction.cpp (modified) (1 diff)
• Actions/MoleculeAction/ForceAnnealingAction.def (modified) (2 diffs)
• Dynamics/ForceAnnealing.hpp (modified) (3 diffs)

src/Actions/MoleculeAction/ForceAnnealingAction.cpp

@@ -117 +117 @@
   // perform optimization step
   LOG(1, "Structural optimization.");
-  optimizer(CurrentStep, 1);
+  optimizer(CurrentStep, 1, params.UseBondGraph.get());
   STATUS("Successfully optimized structure by one step.");
 

src/Actions/MoleculeAction/ForceAnnealingAction.def

@@ -17 +17 @@
 // 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)(int)(double)
-#define paramtokens ("forces-file")("steps")("output-every-step")("max-distance")("damping-factor")
-#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")("damping factor for decreasing the shift with bond graph distance from the causing site")
-#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(PARAM_DEFAULT("0"))(PARAM_DEFAULT(0))(PARAM_DEFAULT(0.5))
-#define paramreferences (forcesfile)(steps)(DoOutput)(MaxDistance)(DampingFactor)
+#define paramtypes (boost::filesystem::path)(unsigned int)(bool)(int)(double)(bool)
+#define paramtokens ("forces-file")("steps")("output-every-step")("max-distance")("damping-factor")("use-bondgraph")
+#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")("damping factor for decreasing the shift with bond graph distance from the causing site")("whether annealing should take place focusing on atoms alone or on the bonds (faster)")
+#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(PARAM_DEFAULT("0"))(PARAM_DEFAULT(0))(PARAM_DEFAULT(0.5))(PARAM_DEFAULT(0))
+#define paramreferences (forcesfile)(steps)(DoOutput)(MaxDistance)(DampingFactor)(UseBondGraph)
 #define paramvalids \
 (DummyValidator< boost::filesystem::path >()) \
@@ -27 +27 @@
 (DummyValidator<bool>()) \
 (DummyValidator< int >()) \
-(RangeValidator< double >(0,1))
+(RangeValidator< double >(0,1)) \
+(DummyValidator<bool>())
 
 #define statetypes (std::vector<AtomicInfo>)(std::vector<AtomicInfo>)

src/Dynamics/ForceAnnealing.hpp

@@ -90 +90 @@
    *
    *
-   * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param CurrentTimeStep 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 size_t offset)
+  void operator()(
+      const int _CurrentTimeStep,
+      const size_t _offset,
+      const bool _UseBondgraph)
   {
     // make sum of forces equal zero
-    AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(offset,CurrentTimeStep);
+    AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(_offset, _CurrentTimeStep);
 
     // are we in initial step? Then set static entities
+    Vector maxComponents(zeroVec);
     if (currentStep == 0) {
       currentDeltat = AtomicForceManipulator<T>::Deltat;
       currentStep = 1;
       LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
+
+      // always use atomic annealing on first step
+      anneal(_CurrentTimeStep, _offset, maxComponents);
     } else {
       ++currentStep;
       LOG(2, "DEBUG: current step is #" << currentStep);
-    }
-
+
+      if (_UseBondgraph)
+        annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
+      else
+        anneal(_CurrentTimeStep, _offset, maxComponents);
+    }
+
+
+    LOG(1, "STATUS: Largest remaining force components at step #"
+        << currentStep << " are " << maxComponents);
+
+    // are we in final step? Remember to reset static entities
+    if (currentStep == maxSteps) {
+      LOG(2, "DEBUG: Final step, resetting values");
+      reset();
+    }
+  }
+
+  /** Performs Gradient optimization on the atoms.
+   *
+   * 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 offset offset in matrix file to the first force component
+   * \param maxComponents to be filled with maximum force component over all atoms
+   */
+  void anneal(
+      const int CurrentTimeStep,
+      const size_t offset,
+      Vector &maxComponents)
+  {
+                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(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
+                        const Vector currentPosition = (*iter)->getPosition();
+                        const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
+                        const Vector currentGradient = (*iter)->getAtomicForce();
+      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);
+//                      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;
+                        }
+      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];
+
+                        // 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(), CurrentTimeStep-1, Deltat, IsAngstroem);
+
+                        // reset force vector for next step except on final one
+                        if (currentStep != maxSteps)
+                          (*iter)->setAtomicForce(zeroVec);
+                }
+
+                LOG(1, "STATUS: Largest remaining force components at step #"
+                    << currentStep << " are " << maxComponents);
+
+                // are we in final step? Remember to reset static entities
+                if (currentStep == maxSteps) {
+                  LOG(2, "DEBUG: Final step, resetting values");
+                  reset();
+                }
+  }
+
+  /** Performs Gradient optimization on the bonds.
+   *
+   * We assume that forces have just been calculated. These forces are projected
+   * onto the bonds and these are annealed subsequently by moving atoms in the
+   * bond neighborhood on either side conjunctively.
+   *
+   *
+   * \param CurrentTimeStep 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
+   * \param maxComponents to be filled with maximum force component over all atoms
+   */
+  void annealWithBondGraph(
+      const int CurrentTimeStep,
+      const size_t offset,
+      Vector &maxComponents)
+  {
     // get nodes on either side of selected bond via BFS discovery
     BoostGraphCreator BGcreator;
@@ -375 +494 @@
     }
 
-    Vector maxComponents(zeroVec);
     for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
         iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
@@ -401 +519 @@
       walker->setAtomicVelocity(update);
       LOG(3, "DEBUG: Applying update " << update << " to atom " << *walker);
-    }
-
-    LOG(1, "STATUS: Largest remaining force components at step #"
-        << currentStep << " are " << maxComponents);
-
-    // are we in final step? Remember to reset static entities
-    if (currentStep == maxSteps) {
-      LOG(2, "DEBUG: Final step, resetting values");
-      reset();
     }
   }