Changeset 6d360a for src/Dynamics/ForceAnnealing.hpp

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.

File:

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

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