Context Navigation

← Previous Changeset
Next Changeset →

Changeset 0682c2

Timestamp:

Jun 20, 2018, 8:21:41 AM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph

Children:

3183f5

Parents:

bd19c1

git-author:

Frederik Heber <frederik.heber@…> (06/16/18 00:03:03)

git-committer:

Frederik Heber <frederik.heber@…> (06/20/18 08:21:41)

Message:

Excluded hydrogens due to their small mass from spreading updates.

hydrogens will rather adapt to position changes of other atoms than the other way round.

File:

: 1 edited

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

Legend:

: Unmodified
: Added
: Removed

src/Dynamics/ForceAnnealing.hpp

-              rbd19c1
+              r0682c2
+  }
+        /** Helper function to insert \a PositionUpdate into a map for every atom.
+         *
+         * \param _GatheredUpdates map of updates per atom
+         * \param _LargestUpdate_per_Atom map with the largest update per atom for checking
+         * \param _atomno key for map
+         * \param _PositionUpdate update to add
+         * \param _factor optional dampening factor
+         */
+        void updateInserter(
+            std::map<atomId_t, Vector> &_GatheredUpdates,
+            std::map<atomId_t, double> &_LargestUpdate_per_Atom,
+            const atomId_t _atomno,
+            const Vector &_PositionUpdate,
+            const double _factor = 1.
+            )
+        {
+          if (_GatheredUpdates.count(_atomno)) {
+            _GatheredUpdates[_atomno] += _factor*_PositionUpdate;
+            _LargestUpdate_per_Atom[_atomno] =
+                std::max(_LargestUpdate_per_Atom[_atomno], _factor*_PositionUpdate.Norm());
+          } else {
+            _GatheredUpdates.insert(
+                std::make_pair(
+                    _atomno,
+                    _factor*_PositionUpdate) );
+            _LargestUpdate_per_Atom.insert(
+                std::make_pair(
+                    _atomno,
+                    _PositionUpdate.Norm()) );
+          }
+        }
   /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
+   *
 …
       LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
+      /** for each atom, we imagine a plane at the position of the atom with
+       * its atomic gradient as the normal vector. We go through all its bonds
+       * and check on which side of the plane the bond is. This defines whether
+       * the bond is contracting (+) or expanding (-) with respect to this atom.
+       *
+       * A bond has two atoms, however. Hence, we do this for either atom and
+       * look at the combination: Is it in sum contracting or expanding given
+       * both projected_forces?
+       */
+      /** go through all bonds and check projected_forces and side of plane
+       * the idea is that if all bonds on one side are contracting ones or expanding,
+       * respectively, then we may shift not only the atom with respect to its
+       * gradient but also its neighbors (towards contraction or towards
+       * expansion depending on direction of gradient).
+       * if they are mixed on both sides of the plane, then we simply shift
+       * only the atom itself.
+       * if they are not mixed on either side, then we also only shift the
+       * atom, namely away from expanding and towards contracting bonds.
+       *
+       * We may get this information right away by looking at the projected_forces.
+       * They give the atomic gradient of either atom projected onto the BondVector
+       * with an additional weight in [0,1].
+       */
+      // sign encodes side of plane and also encodes contracting(-) or expanding(+)
+      typedef std::vector<int> sides_t;
+      typedef std::vector<int> types_t;
+      sides_t sides;
+      types_t types;
+      const BondList& ListOfBonds = walker.getListOfBonds();
+      for(BondList::const_iterator bonditer = ListOfBonds.begin();
+          bonditer != ListOfBonds.end(); ++bonditer) {
+        const bond::ptr &current_bond = *bonditer;
+        // BondVector goes from bond::rightatom to bond::leftatom
+        const size_t index = bv.getIndexForBond(current_bond);
+        std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
+            projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside];
+        // note that projected_forces has sign such as to indicate whether
+        // atomic gradient wants bond to contract (-) or expand (+).
+        // This goes into sides: Minus side points away from gradient, plus side point
+        // towards gradient.
+        //
+        // the sum of both bond sides goes into types, depending on which is
+        // stronger if either wants a different thing
+        const double &temp = forcelist[index];
+        if (fabs(temp) < MYEPSILON)
+          sides.push_back(1);
+        else
+          sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision
+        ASSERT( (sides.back() == 1) || (sides.back() == -1),
+            "ForceAnnealing() - sides is not in {-1,1}.");
+        const double sum =
+            projected_forces[BondVectors::leftside][index]+projected_forces[BondVectors::rightside][index];
+        types.push_back( sum/fabs(sum) );
+        LOG(4, "DEBUG: Bond " << *current_bond << " is on side " << sides.back()
+            << " and has type " << types.back());
+      }
+//      /// check whether both conditions are compatible:
+//      // i.e. either we have ++/-- for all entries in sides and types
+//      // or we have +-/-+ for all entries
+//      // hence, multiplying and taking the sum and its absolute value
+//      // should be equal to the maximum number of entries
+//      sides_t results;
+//      std::transform(
+//          sides.begin(), sides.end(),
+//          types.begin(),
+//          std::back_inserter(results),
+//          std::multiplies<int>);
+//      int result = abs(std::accumulate(results.begin(), results.end(), 0, std::plus<int>));
+      std::vector<size_t> first_per_side(2, (size_t)-1); //!< mark down one representative from either side
+      std::vector< std::vector<int> > types_per_side(2); //!< gather all types on each side
+      types_t::const_iterator typesiter = types.begin();
+      for (sides_t::const_iterator sidesiter = sides.begin();
+          sidesiter != sides.end(); ++sidesiter, ++typesiter) {
+        const size_t index = (*sidesiter+1)/2;
+        types_per_side[index].push_back(*typesiter);
+        if (first_per_side[index] == (size_t)-1)
+          first_per_side[index] = std::distance(const_cast<const sides_t &>(sides).begin(), sidesiter);
+      }
+      LOG(4, "DEBUG: First on side minus is " << first_per_side[0] << ", and first on side plus is "
+          << first_per_side[1]);
+      //!> enumerate types per side with a little witching with the numbers to allow easy setting from types
+      enum whichtypes_t {
+        contracting=0,
+        unset=1,
+        expanding=2,
+        mixed
+      };
+      std::vector<int> typeside(2, unset);
+      for(size_t i=0;i<2;++i) {
+        for (std::vector<int>::const_iterator tpsiter = types_per_side[i].begin();
+            tpsiter != types_per_side[i].end(); ++tpsiter) {
+          if (typeside[i] == unset) {
+            typeside[i] = *tpsiter+1; //contracting(0) or expanding(2)
+          } else {
+            if (typeside[i] != (*tpsiter+1)) // no longer he same type
+              typeside[i] = mixed;
+      if (walker.getElementNo() != 1) {
+        /** for each atom, we imagine a plane at the position of the atom with
+         * its atomic gradient as the normal vector. We go through all its bonds
+         * and check on which side of the plane the bond is. This defines whether
+         * the bond is contracting (+) or expanding (-) with respect to this atom.
+         *
+         * A bond has two atoms, however. Hence, we do this for either atom and
+         * look at the combination: Is it in sum contracting or expanding given
+         * both projected_forces?
+         */
+        /** go through all bonds and check projected_forces and side of plane
+         * the idea is that if all bonds on one side are contracting ones or expanding,
+         * respectively, then we may shift not only the atom with respect to its
+         * gradient but also its neighbors (towards contraction or towards
+         * expansion depending on direction of gradient).
+         * if they are mixed on both sides of the plane, then we simply shift
+         * only the atom itself.
+         * if they are not mixed on either side, then we also only shift the
+         * atom, namely away from expanding and towards contracting bonds.
+         *
+         * We may get this information right away by looking at the projected_forces.
+         * They give the atomic gradient of either atom projected onto the BondVector
+         * with an additional weight in [0,1].
+         */
+        // sign encodes side of plane and also encodes contracting(-) or expanding(+)
+        typedef std::vector<int> sides_t;
+        typedef std::vector<int> types_t;
+        sides_t sides;
+        types_t types;
+        const BondList& ListOfBonds = walker.getListOfBonds();
+        for(BondList::const_iterator bonditer = ListOfBonds.begin();
+            bonditer != ListOfBonds.end(); ++bonditer) {
+          const bond::ptr &current_bond = *bonditer;
+          // BondVector goes from bond::rightatom to bond::leftatom
+          const size_t index = bv.getIndexForBond(current_bond);
+          std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
+              projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside];
+          // note that projected_forces has sign such as to indicate whether
+          // atomic gradient wants bond to contract (-) or expand (+).
+          // This goes into sides: Minus side points away from gradient, plus side point
+          // towards gradient.
+          //
+          // the sum of both bond sides goes into types, depending on which is
+          // stronger if either wants a different thing
+          const double &temp = forcelist[index];
+          if (fabs(temp) < MYEPSILON)
+            sides.push_back(1);
+          else
+            sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision
+          ASSERT( (sides.back() == 1) || (sides.back() == -1),
+              "ForceAnnealing() - sides is not in {-1,1}.");
+          const double sum =
+              projected_forces[BondVectors::leftside][index]+projected_forces[BondVectors::rightside][index];
+          types.push_back( sum/fabs(sum) );
+          LOG(4, "DEBUG: Bond " << *current_bond << " is on side " << sides.back()
+              << " and has type " << types.back());
+        }
+  //      /// check whether both conditions are compatible:
+  //      // i.e. either we have ++/-- for all entries in sides and types
+  //      // or we have +-/-+ for all entries
+  //      // hence, multiplying and taking the sum and its absolute value
+  //      // should be equal to the maximum number of entries
+  //      sides_t results;
+  //      std::transform(
+  //          sides.begin(), sides.end(),
+  //          types.begin(),
+  //          std::back_inserter(results),
+  //          std::multiplies<int>);
+  //      int result = abs(std::accumulate(results.begin(), results.end(), 0, std::plus<int>));
+        std::vector<size_t> first_per_side(2, (size_t)-1); //!< mark down one representative from either side
+        std::vector< std::vector<int> > types_per_side(2); //!< gather all types on each side
+        types_t::const_iterator typesiter = types.begin();
+        for (sides_t::const_iterator sidesiter = sides.begin();
+            sidesiter != sides.end(); ++sidesiter, ++typesiter) {
+          const size_t index = (*sidesiter+1)/2;
+          types_per_side[index].push_back(*typesiter);
+          if (first_per_side[index] == (size_t)-1)
+            first_per_side[index] = std::distance(const_cast<const sides_t &>(sides).begin(), sidesiter);
+        }
+        LOG(4, "DEBUG: First on side minus is " << first_per_side[0] << ", and first on side plus is "
+            << first_per_side[1]);
+        //!> enumerate types per side with a little witching with the numbers to allow easy setting from types
+        enum whichtypes_t {
+          contracting=0,
+          unset=1,
+          expanding=2,
+          mixed
+        };
+        std::vector<int> typeside(2, unset);
+        for(size_t i=0;i<2;++i) {
+          for (std::vector<int>::const_iterator tpsiter = types_per_side[i].begin();
+              tpsiter != types_per_side[i].end(); ++tpsiter) {
+            if (typeside[i] == unset) {
+              typeside[i] = *tpsiter+1; //contracting(0) or expanding(2)
+            } else {
+              if (typeside[i] != (*tpsiter+1)) // no longer he same type
+                typeside[i] = mixed;
+            }
+          }
+        }
+      }
+      LOG(4, "DEBUG: Minus side is " << typeside[0] << " and plus side is " << typeside[1]);
       typedef std::vector< std::pair<atomId_t, atomId_t> > RemovedEdges_t;
       if ((typeside[0] != mixed) || (typeside[1] != mixed)) {
         const size_t sideno = ((typeside[0] != mixed) && (typeside[0] != unset)) ? 0 : 1;
         LOG(4, "DEBUG: Chosen side is " << sideno << " with type " << typeside[sideno]);
         ASSERT( (typeside[sideno] == contracting) || (typeside[sideno] == expanding),
             "annealWithBondGraph_BB() - chosen side is neither expanding nor contracting.");
         // one side is not mixed, all bonds on one side are of same type
         // hence, find out which bonds to exclude
+        const BondList& ListOfBonds = walker.getListOfBonds();
         // sideno is away (0) or in direction (1) of gradient
         // tpyes[first_per_side[sideno]] is either contracting (-1) or expanding (+1)
         // : side (i), where (i) means which bonds we keep for the BFS, bonds
         // on side (-i) are removed
         // If all bonds on side away (0) want expansion (+1), move towards side with atom: side 1
+        // if all bonds side towards (1) want contraction (-1), move away side with atom : side -1
         // unsure whether this or do nothing in the remaining cases:
         // If all bonds on side toward (1) want expansion (+1), move away side with atom : side -1
         //    (the reasoning is that the bond's other atom must have a stronger
         //     gradient in the same direction and they push along atoms in
         //     gradient direction: we don't want to interface with those.
         //     Hence, move atoms along on away side
         // if all bonds side away (0) want contraction (-1), move towards side with atom: side 1
         //    (the reasoning is the same, don't interfere with update from
         //     stronger gradient)
         // hence, the decision is only based on sides once we have picked a side
+        // depending on all bonds associated with have same good type.
         // away from gradient (minus) and contracting
         // or towards gradient (plus) and expanding
         // gather all on same side and remove
         const double sign =
             (sides[first_per_side[sideno]] == types[first_per_side[sideno]])
+            ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]];
         LOG(4, "DEBUG: Removing edges from side with sign " << sign);
         BondList::const_iterator bonditer = ListOfBonds.begin();
         RemovedEdges_t RemovedEdges;
         for (sides_t::const_iterator sidesiter = sides.begin();
             sidesiter != sides.end(); ++sidesiter, ++bonditer) {
           if (*sidesiter == sign) {
             // remove the edge
             const bond::ptr &current_bond = *bonditer;
             LOG(5, "DEBUG: Removing edge " << *current_bond);
             RemovedEdges.push_back( std::make_pair(
                 current_bond->leftatom->getId(),
                 current_bond->rightatom->getId())
             );
+#ifndef NDEBUG
             const bool status =
+#endif
                 BGcreator.removeEdge(RemovedEdges.back());
             ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
+        LOG(4, "DEBUG: Minus side is " << typeside[0] << " and plus side is " << typeside[1]);
+        typedef std::vector< std::pair<atomId_t, atomId_t> > RemovedEdges_t;
+        if ((typeside[0] != mixed) || (typeside[1] != mixed)) {
+          const size_t sideno = ((typeside[0] != mixed) && (typeside[0] != unset)) ? 0 : 1;
+          LOG(4, "DEBUG: Chosen side is " << sideno << " with type " << typeside[sideno]);
+          ASSERT( (typeside[sideno] == contracting) || (typeside[sideno] == expanding),
+              "annealWithBondGraph_BB() - chosen side is neither expanding nor contracting.");
+          // one side is not mixed, all bonds on one side are of same type
+          // hence, find out which bonds to exclude
+          const BondList& ListOfBonds = walker.getListOfBonds();
+          // sideno is away (0) or in direction (1) of gradient
+          // tpyes[first_per_side[sideno]] is either contracting (-1) or expanding (+1)
+          // : side (i), where (i) means which bonds we keep for the BFS, bonds
+          // on side (-i) are removed
+          // If all bonds on side away (0) want expansion (+1), move towards side with atom: side 1
+          // if all bonds side towards (1) want contraction (-1), move away side with atom : side -1
+          // unsure whether this or do nothing in the remaining cases:
+          // If all bonds on side toward (1) want expansion (+1), move away side with atom : side -1
+          //    (the reasoning is that the bond's other atom must have a stronger
+          //     gradient in the same direction and they push along atoms in
+          //     gradient direction: we don't want to interface with those.
+          //     Hence, move atoms along on away side
+          // if all bonds side away (0) want contraction (-1), move towards side with atom: side 1
+          //    (the reasoning is the same, don't interfere with update from
+          //     stronger gradient)
+          // hence, the decision is only based on sides once we have picked a side
+          // depending on all bonds associated with have same good type.
+          // away from gradient (minus) and contracting
+          // or towards gradient (plus) and expanding
+          // gather all on same side and remove
+          const double sign =
+              (sides[first_per_side[sideno]] == types[first_per_side[sideno]])
+              ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]];
+          LOG(4, "DEBUG: Removing edges from side with sign " << sign);
+          BondList::const_iterator bonditer = ListOfBonds.begin();
+          RemovedEdges_t RemovedEdges;
+          for (sides_t::const_iterator sidesiter = sides.begin();
+              sidesiter != sides.end(); ++sidesiter, ++bonditer) {
+            if (*sidesiter == sign) {
+              // remove the edge
+              const bond::ptr &current_bond = *bonditer;
+              LOG(5, "DEBUG: Removing edge " << *current_bond);
+              RemovedEdges.push_back( std::make_pair(
+                  current_bond->leftatom->getId(),
+                  current_bond->rightatom->getId())
+              );
+  #ifndef NDEBUG
+              const bool status =
+  #endif
+                  BGcreator.removeEdge(RemovedEdges.back());
+              ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
+            }
+          }
+        }
+        // perform limited-horizon BFS
+        BoostGraphHelpers::Nodeset_t bondside_set;
+        BreadthFirstSearchGatherer::distance_map_t distance_map;
+        bondside_set = NodeGatherer(walker.getId(), max_distance);
+        distance_map = NodeGatherer.getDistances();
+        std::sort(bondside_set.begin(), bondside_set.end());
+        // re-add edge
+        for (RemovedEdges_t::const_iterator edgeiter = RemovedEdges.begin();
+            edgeiter != RemovedEdges.end(); ++edgeiter)
+          BGcreator.addEdge(edgeiter->first, edgeiter->second);
+        // update position with dampening factor on the discovered bonds
+        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
+            setiter != bondside_set.end(); ++setiter) {
+          const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
+            = distance_map.find(*setiter);
+          ASSERT( diter != distance_map.end(),
+              "ForceAnnealing() - could not find distance to an atom.");
+          const double factor = pow(damping_factor, diter->second+1);
+          LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
+              << factor << "*" << PositionUpdate);
+          if (GatheredUpdates.count((*setiter))) {
+            GatheredUpdates[(*setiter)] += factor*PositionUpdate;
+            LargestUpdate_per_Atom[(*setiter)] =
+                std::max(LargestUpdate_per_Atom[(*setiter)], factor*PositionUpdate.Norm());
+          } else {
+            GatheredUpdates.insert(
+                std::make_pair(
+                    (*setiter),
+                    factor*PositionUpdate) );
+            LargestUpdate_per_Atom.insert(
+                std::make_pair(
+                    (*setiter),
+                    factor*PositionUpdate.Norm()) );
+          // perform limited-horizon BFS
+          BoostGraphHelpers::Nodeset_t bondside_set;
+          BreadthFirstSearchGatherer::distance_map_t distance_map;
+          bondside_set = NodeGatherer(walker.getId(), max_distance);
+          distance_map = NodeGatherer.getDistances();
+          std::sort(bondside_set.begin(), bondside_set.end());
+          // re-add edge
+          for (RemovedEdges_t::const_iterator edgeiter = RemovedEdges.begin();
+              edgeiter != RemovedEdges.end(); ++edgeiter)
+            BGcreator.addEdge(edgeiter->first, edgeiter->second);
+          // update position with dampening factor on the discovered bonds
+          for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
+              setiter != bondside_set.end(); ++setiter) {
+            const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
+              = distance_map.find(*setiter);
+            ASSERT( diter != distance_map.end(),
+                "ForceAnnealing() - could not find distance to an atom.");
+            const double factor = pow(damping_factor, diter->second+1);
+            LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
+                << factor << "*" << PositionUpdate);
+            updateInserter(GatheredUpdates, LargestUpdate_per_Atom, *setiter, PositionUpdate, factor);
+          }
+        } else {
+          // simple atomic annealing, i.e. damping factor of 1
+          updateInserter(GatheredUpdates, LargestUpdate_per_Atom, walker.getId(), PositionUpdate);
+        }
       } else {
+        // simple atomic annealing, i.e. damping factor of 1
+        LOG(3, "DEBUG: Update for atom #" << walker.getId() << " will be " << PositionUpdate);
+        GatheredUpdates.insert(
+            std::make_pair(
+                walker.getId(),
+                PositionUpdate) );
+        LargestUpdate_per_Atom.insert(
+            std::make_pair(
+                walker.getId(),
+                PositionUpdate.Norm()) );
+        // hydrogens (are light-weighted and therefore) are always updated normally
+        LOG(3, "DEBUG: Update for hydrogen #" << walker.getId() << " will be " << PositionUpdate);
+        updateInserter(GatheredUpdates, LargestUpdate_per_Atom, walker.getId(), PositionUpdate);
+      }
+    }

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 0682c2

Legend:

src/Dynamics/ForceAnnealing.hpp

Download in other formats: