Changeset 698308 for src/Dynamics

Timestamp:

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

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_with_BondGraph_continued

Children:

12d20c

Parents:

b7fe6e

git-author:

Frederik Heber <frederik.heber@…> (06/27/17 21:17:18)

git-committer:

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

Message:

Fixing annealWithBondgraph().

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -441 +441 @@
         const double &bondforce = projected_forces[0][side][index];
         const double &oldbondforce = projected_forces[1][side][index];
-        const double bondforcedifference = (bondforce - oldbondforce);
+        const double bondforcedifference = fabs(bondforce - oldbondforce);
+        LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
+            << " side of bond is " << bondforce);
+        LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
+            << " side of bond is " << oldbondforce);
         // if difference or bondforce itself is zero, do nothing
         if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
           continue;
+
+        // get BondVector to bond
         const BondVectors::mapped_t::const_iterator bviter =
             bondvectors.find(current_bond);
         ASSERT( bviter != bondvectors.end(),
             "ForceAnnealing() - cannot find current_bond ?");
+        ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
+            "ForceAnnealing() - norm of BondVector is not one");
         const Vector &BondVector = bviter->second;
-        const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
-        const Vector &currentPosition = bondatom[side]->getPosition();
-        const double stepwidth =
-            fabs((currentPosition - oldPosition).ScalarProduct(BondVector))/bondforcedifference;
-        Vector PositionUpdate = stepwidth * BondVector;
+
+        // calculate gradient and position differences for stepwidth
+        const Vector currentGradient = bondforce * BondVector;
+        LOG(4, "DEBUG: current projected gradient for "
+            << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
+        const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
+        const Vector &currentPosition = bondatom[side]->getPositionAtStep(currentStep);
+        const Vector PositionDifference = currentPosition - oldPosition;
+        LOG(4, "DEBUG: old position is " << oldPosition);
+        LOG(4, "DEBUG: current position is " << currentPosition);
+        LOG(4, "DEBUG: difference in position is " << PositionDifference);
+        LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
+        LOG(4, "DEBUG: bondvector is " << BondVector);
+
+        // calculate step width
+        double stepwidth =
+            fabs(PositionDifference.ScalarProduct(BondVector))/bondforcedifference;
         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 * BondVector;
+          stepwidth = currentDeltat;
         }
-        LOG(3, "DEBUG: Update would be " << PositionUpdate);
+        Vector PositionUpdate = stepwidth * currentGradient;
+        LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
 
         // add PositionUpdate for all nodes in the bondside_set
@@ -470 +491 @@
           ASSERT( diter != distance_map[side].end(),
               "ForceAnnealing() - could not find distance to an atom.");
-          const double factor = pow(damping_factor, diter->second);
+          const double factor = pow(damping_factor, diter->second+1);
           LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
               << factor << "*" << PositionUpdate);