Changeset bd0159 for src/Dynamics/ForceAnnealing.hpp

Timestamp:

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

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_oldresults

Children:

8754b1

Parents:

84f449

git-author:

Frederik Heber <heber@…> (04/06/17 05:09:37)

git-committer:

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

Message:

ForceAnnealing now uses step width according to Barzilai-Borwein method.

• this is not as good as what MPQC does but as long as Actions don't have an internal state, we cannot do anything better, e.g. a real line-search.
• removed paramater delta from optimize-structure.
• DOCU: Updated entry in user guide and noted Barzilai-Borwein method.
• TESTS: removed deltat from force annealing regression tests.

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -32 +32 @@
  *
  * 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
+ * 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
  * calculation of the forces).
@@ -42 +42 @@
   /** Constructor of class ForceAnnealing.
    *
+   * \note We use a fixed delta t of 1.
+   *
    * \param _atoms set of atoms to integrate
    * \param _Deltat time step width in atomic units
@@ -49 +51 @@
   ForceAnnealing(
       AtomSetMixin<T> &_atoms,
-      double _Deltat,
       bool _IsAngstroem,
       const size_t _maxSteps) :
-    AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
+    AtomicForceManipulator<T>(_atoms, 1., _IsAngstroem),
     maxSteps(_maxSteps)
   {}
@@ -89 +90 @@
         iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
       // atom's force vector gives steepest descent direction
+      const Vector oldPosition = (*iter)->getPositionAtStep(NextStep-2 >= 0 ? NextStep - 2 : 0);
       const Vector currentPosition = (*iter)->getPosition();
+      const Vector oldGradient = (*iter)->getAtomicForceAtStep(NextStep-2 >= 0 ? NextStep - 2 : 0);
       const Vector currentGradient = (*iter)->getAtomicForce();
       LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
 
-      // artificial update: deltat may be considered as 1/2 s^2 units, mass
-      // is neglected deliberately as this makes all atoms equally fast or
-      // hydrogens slower (and they need to wait for other atoms to arrive at
-      // final position).
-      Vector PositionUpdate = currentDeltat * currentGradient;
+      // we use Barzilai-Borwein update with position reversed to get descent
+      const Vector GradientDifference = (currentGradient - oldGradient);
+      const double stepwidth =
+          fabs((currentPosition - oldPosition).ScalarProduct(GradientDifference))/
+          GradientDifference.NormSquared();
+      Vector PositionUpdate = stepwidth * currentGradient;
+      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.");
+        PositionUpdate = currentDeltat * currentGradient;
+      }
       LOG(3, "DEBUG: Update would be " << PositionUpdate);