Changeset efd020 for src

Timestamp:

Apr 10, 2018, 6:43:12 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:

8c6b68

Parents:

12f16c

git-author:

Frederik Heber <frederik.heber@…> (09/11/17 17:03:27)

git-committer:

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

Message:

Readded step-world-time to StructuralOptimization.

• adapted ForceAnnealing and its Action accordingly, i.e. they receive the old time step as the current step.
• TESTFIX: needed to add "--step-world-time 1" to ForceAnnealing regression test, results then remained exactly the same.

Location:

src

Files:

• Actions/FragmentationAction/StructuralOptimizationAction.cpp (modified) (1 diff)
• Actions/MoleculeAction/ForceAnnealingAction.cpp (modified) (2 diffs)
• Dynamics/ForceAnnealing.hpp (modified) (4 diffs)

src/Actions/FragmentationAction/StructuralOptimizationAction.cpp

@@ -69 +69 @@
   prototype_actions.addAction(AR.getActionByName(std::string("fragment-automation")).clone());
   prototype_actions.addAction(AR.getActionByName(std::string("analyse-fragment-results")).clone());
+  prototype_actions.addAction(AR.getActionByName(std::string("step-world-time")).clone());
   prototype_actions.addAction(AR.getActionByName(std::string("force-annealing")).clone());
   prototype_actions.addAction(AR.getActionByName(std::string("output")).clone());

src/Actions/MoleculeAction/ForceAnnealingAction.cpp

@@ -94 +94 @@
       params.DampingFactor.get());
   size_t CurrentStep = WorldTime::getInstance().getTime();
+  if (CurrentStep == 0) {
+    ELOG(1, "WorldTime must be at least at step 1 already, use step-world-time if necessary.");
+    return Action::failure;
+  }
 
-  // parse forces into current step
+  // parse forces into last step (assuming we stepped on already)
   if (!params.forcesfile.get().string().empty()) {
     LOG(1, "Parsing forces file.");
-    if (!optimizer.parseForcesFile(params.forcesfile.get().string().c_str(), CurrentStep))
+    if (!optimizer.parseForcesFile(params.forcesfile.get().string().c_str(), CurrentStep-1))
       LOG(2, "File " << params.forcesfile.get() << " not found.");
     else
@@ -104 +108 @@
   }
 
-  // copy current time step to new one and and proceed on this one
-  {
-    for (World::AtomSelectionIterator iter = World::getInstance().beginAtomSelection();
-        iter != World::getInstance().endAtomSelection();
-        ++iter) {
-      atom * const Walker = iter->second;
-      Walker->setPositionAtStep(CurrentStep+1,
-          Walker->getPositionAtStep(CurrentStep));
-      // force have already been calculated, hence copy them
-      Walker->setAtomicVelocityAtStep(CurrentStep+1,
-          Walker->getAtomicVelocityAtStep(CurrentStep));
-      Walker->setAtomicForceAtStep(CurrentStep+1,
-          Walker->getAtomicForceAtStep(CurrentStep));
-    }
-    // increment to next time step: re-creates bond graph
-    ++CurrentStep;
-    World::getInstance().setTime(CurrentStep);
-  }
-
   // perform optimization step
   LOG(1, "Structural optimization.");
-  optimizer(CurrentStep, 1, params.UseBondGraph.get());
+  optimizer(CurrentStep-1, 1, params.UseBondGraph.get());
   STATUS("Successfully optimized structure by one step.");
 

src/Dynamics/ForceAnnealing.hpp

@@ -139 +139 @@
         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();
+      const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
+      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
+      const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
+      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
       LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
       LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
@@ -197 +197 @@
    *
    *
-   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param CurrentTimeStep current time step (i.e. t where \f$ t + \Delta t \f$ is 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
@@ -226 +226 @@
         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();
+      const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
+      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
+      const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
+      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
       LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
 
@@ -269 +269 @@
         const bond &current_bond = *(*bonditer);
         LOG(2, "DEBUG: Looking at bond " << current_bond);
-        Vector BondVector = (*iter)->getPosition();
+        Vector BondVector = (*iter)->getPositionAtStep(CurrentTimeStep);
         BondVector -= ((*iter)->getId() == current_bond.rightatom->getId())
-            ? current_bond.rightatom->getPosition() : current_bond.leftatom->getPosition();
+            ? current_bond.rightatom->getPositionAtStep(CurrentTimeStep) : current_bond.leftatom->getPositionAtStep(CurrentTimeStep);
         BondVector.Normalize();
         if (BondVector.ScalarProduct(currentGradient) < 0) {