Changeset b2acca for src/Dynamics

Timestamp:

Apr 10, 2018, 6:43:11 AM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Enhanced_StructuralOptimization_continued, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_contraction-expansion, Gui_displays_atomic_force_velocity, PythonUI_with_named_parameters, StoppableMakroAction, TremoloParser_IncreasedPrecision

Children:

216840

Parents:

ecfcf6

git-author:

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

git-committer:

Frederik Heber <frederik.heber@…> (04/10/18 06:43:11)

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.
• TESTS: Removed XFAIL from ForceAnnealing regression tests as they work now again.

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -80 +80 @@
    *
    *
-   * \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 NextStep, 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,NextStep);
+    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
 //    std::vector<atomId_t> atomids;
@@ -116 +235 @@
 
     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();
         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 oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
       const Vector currentPosition = (*iter)->getPosition();
-      const Vector oldGradient = (*iter)->getAtomicForceAtStep(NextStep-2 >= 0 ? NextStep - 2 : 0);
+      const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
       const Vector currentGradient = (*iter)->getAtomicForce();
       LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
@@ -225 +343 @@
       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();
-    }
   }