Changeset bccbe9 for src/Actions

Timestamp:

Jul 12, 2017, 7:10:31 PM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

Action_Thermostats, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Exclude_Hydrogens_annealWithBondGraph, Fix_Verbose_Codepatterns, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, Gui_displays_atomic_force_velocity, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, PythonUI_with_named_parameters, Recreated_GuiChecks, StoppableMakroAction, TremoloParser_IncreasedPrecision

Children:

6e5b8d

Parents:

d24ef58

git-author:

Frederik Heber <frederik.heber@…> (05/18/17 17:45:47)

git-committer:

Frederik Heber <frederik.heber@…> (07/12/17 19:10:31)

Message:

Extracted extraction (subset of) nodes from BoostGraph into BreadthFirstSearchGatherer.

• also added helper namespace BoostGraphHelpers.
• we now treat the vertex indices and vertex names properly. Before that they had to coincide. Now, the name is the atomic id associated with the node and the index is the boost::graph internal index.

Location:

src/Actions/MoleculeAction

Files:

• StretchBondAction.cpp (modified) (8 diffs)
• StretchBondAction.def (modified) (1 diff)

src/Actions/MoleculeAction/StretchBondAction.cpp

@@ -33 +33 @@
 #endif
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/breadth_first_search.hpp>
-
 //#include "CodePatterns/MemDebug.hpp"
 
@@ -52 +49 @@
 #include "Descriptors/AtomIdDescriptor.hpp"
 #include "Graph/BoostGraphCreator.hpp"
+#include "Graph/BoostGraphHelpers.hpp"
+#include "Graph/BreadthFirstSearchGatherer.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
@@ -61 +60 @@
 #include "Action_impl_pre.hpp"
 
-
-/**
- * I have no idea why this is so complicated with BGL ...
- *
- * This is taken from the book "The Boost Graph Library: User Guide and Reference Manual, Portable Documents",
- * chapter "Basic Graph Algorithms", example on calculating the bacon number.
- */
-template <typename DistanceMap>
-class distance_recorder : public boost::default_bfs_visitor
-{
-public:
-  distance_recorder(DistanceMap dist) : d(dist) {}
-
-  template <typename Edge, typename Graph>
-  void tree_edge(Edge e, const Graph &g) const {
-    typename boost::graph_traits<Graph>::vertex_descriptor u = source(e,g), v = target(e,g);
-    d[v] = d[u] + 1;
-  }
-
-private:
-  DistanceMap d;
-};
-
-template <typename DistanceMap>
-distance_recorder<DistanceMap> record_distance(DistanceMap d)
-{
-  return distance_recorder<DistanceMap>(d);
-}
 
 static bool addEdgePredicate(
@@ -110 +81 @@
     return Action::failure;
   }
+  molecule *mol = World::getInstance().
+      getMolecule(MoleculeById(atoms[0]->getMolecule()->getId()));
+  if (mol != atoms[1]->getMolecule()) {
+    STATUS("The two selected atoms must belong to the same molecule.");
+    return Action::failure;
+  }
+
+  // gather undo information
+  const double olddistance = atoms[0]->getPosition().distance(atoms[1]->getPosition());
+  const double newdistance = params.bonddistance.get();
+  LOG(1, "INFO: Old bond distance is " << olddistance << ", stretching to " << newdistance << ".");
+
+  // gather sorted ids
   std::vector<atomId_t> atomids(2);
   atomids[0] = atoms[0]->getId();
@@ -115 +99 @@
   std::sort(atomids.begin(), atomids.end());
   LOG(1, "DEBUG: Selected nodes are " << atomids);
-  molecule *mol = World::getInstance().
-      getMolecule(MoleculeById(atoms[0]->getMolecule()->getId()));
-  if (mol != atoms[1]->getMolecule()) {
-    STATUS("The two selected atoms must belong to the same molecule.");
-    return Action::failure;
-  }
-  // gather undo information
-  const double olddistance = atoms[0]->getPosition().distance(atoms[1]->getPosition());
-  const double newdistance = params.bonddistance.get();
-  LOG(1, "INFO: Old bond distance is " << olddistance << ", stretching to " << newdistance << ".");
+
+  const Vector NormalVector = (atoms[0]->getPosition() - atoms[1]->getPosition())* (1./olddistance);
+  const double shift = 0.5*(newdistance - olddistance);
+  std::vector<Vector> Shift(2);
+  Shift[0] = shift * NormalVector;
+  Shift[1] = -shift * NormalVector;
+  Box &domain = World::getInstance().getDomain();
 
   // Assume the selected bond splits the molecule into two parts, each one on
@@ -133 +114 @@
   // leaving the other positions untouched.
 
-  // convert BondGraph into boost::graph
+  // get nodes on either side of selected bond via BFS discovery
   BoostGraphCreator BGcreator;
   BGcreator.createFromMolecule(*mol,
       boost::bind(addEdgePredicate, _1, boost::ref(atomids)));
-  const BoostGraphCreator::UndirectedGraph &molgraph = BGcreator.get();
-
-  const size_t num_vertices = BGcreator.getNumVertices();
-  std::vector< std::vector<size_t> > distances;
-  for (size_t i=0;i<2;++i) {
-    distances.push_back(std::vector<size_t>(num_vertices, num_vertices+1)); // set distance to num+1
-    distances[i][atomids[i]] = 0;
-    boost::breadth_first_search(
-        molgraph,
-        boost::vertex(atomids[i], molgraph),
-        boost::visitor(record_distance(&(distances[i][0]))));
-    LOG(3, "DEBUG: From atom #" << atomids[i]
-        << " BFS discovered the following distances " << distances[i]);
-  }
-
-  const Vector NormalVector = (atoms[0]->getPosition() - atoms[1]->getPosition())* (1./olddistance);
-  const double shift = 0.5*(newdistance - olddistance);
-  std::vector<Vector> Shift(2);
-  Shift[0] = shift * NormalVector;
-  Shift[1] = -shift * NormalVector;
-  Box &domain = World::getInstance().getDomain();
-  std::vector< std::vector<size_t> > bondside_sets(2);
-
-  // Check whether there are common nodes in each set of distances
-  for (size_t i=0;i<num_vertices;++i)
-    if ((distances[0][i] != (num_vertices+1)) && (distances[1][i] != (num_vertices+1))) {
-      ELOG(2, "Node #" << i << " is reachable from either side of bond, hence must by cyclic."
+  BreadthFirstSearchGatherer NodeGatherer(BGcreator);
+  std::vector< BoostGraphHelpers::Nodeset_t > bondside_sets(2);
+  for(size_t j=0;j<2;++j) {
+    bondside_sets[j] = NodeGatherer(atoms[j]->getId());
+    std::sort(bondside_sets[j].begin(), bondside_sets[j].end());
+  }
+
+  // simple test whether bond has split the system in two disjoint sets or not
+  bool isCyclic = false;
+  if ((bondside_sets[0].size() + bondside_sets[1].size()) > BGcreator.getNumVertices()) {
+      // Check whether there are common nodes in each set of distances
+    if (BoostGraphHelpers::isCommonNodeInVector(bondside_sets[0], bondside_sets[1])) {
+      ELOG(2, "Sets contain common node, hence bond must have been by cyclic."
           << " Shifting only bond partners.");
       for(size_t j=0;j<2;++j) {
-        bondside_sets[j].push_back(atoms[j]->getId());
+        bondside_sets[j].clear();
+        bondside_sets[j].push_back( atomids[j] );
         const Vector &position = atoms[j]->getPosition();
         atoms[j]->setPosition( domain.enforceBoundaryConditions(position+Shift[j]) );
       }
-      break;
-    }
+      isCyclic = true;
+    }
+  }
 
   // go through the molecule and stretch each atom in either set of nodes
-  if (bondside_sets[0].empty()) {
-    const BoostGraphCreator::index_map_t index_map = BGcreator.getIndexMap();
+  if (!isCyclic) {
     for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
       const Vector &position = (*iter)->getPosition();
       // for each atom determine in which set of nodes it is and shift accordingly
-      size_t i=0;
-      const size_t nodeindex = boost::get(index_map, boost::vertex((*iter)->getId(), molgraph));
-      for (;i<2;++i) {
-        if (distances[i][nodeindex] != (num_vertices+1)) {
-          (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[i]) );
-          bondside_sets[i].push_back((*iter)->getId());
-          break;
-        }
-      }
-      if (i==2) {
+      const atomId_t &atomid = (*iter)->getId();
+      if (std::binary_search(bondside_sets[0].begin(), bondside_sets[0].end(), atomid)) {
+        (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[0]) );
+      } else if (std::binary_search(bondside_sets[1].begin(), bondside_sets[1].end(), atomid)) {
+        (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[1]) );
+      } else {
         ELOG(1, "Atom " << *iter << " is not contained on either side of bond? Undoing done shifts");
         // Have to undo shifts
-        for (i=0;i<2;++i) {
-          for (std::vector<size_t>::const_iterator iter = bondside_sets[i].begin();
+        for (size_t i=0;i<2;++i) {
+          for (BoostGraphHelpers::Nodeset_t::const_iterator iter = bondside_sets[i].begin();
               iter != bondside_sets[i].end(); ++iter) {
             atom &walker = *World::getInstance().getAtom(AtomById(*iter));
@@ -214 +178 @@
   Box &domain = World::getInstance().getDomain();
   for (size_t i=0;i<2;++i) {
-    for (std::vector<size_t>::const_iterator iter = state->bondside_sets[i].begin();
+    for (BoostGraphHelpers::Nodeset_t::const_iterator iter = state->bondside_sets[i].begin();
         iter != state->bondside_sets[i].end(); ++iter) {
       atom &walker = *World::getInstance().getAtom(AtomById(*iter));
@@ -230 +194 @@
   Box &domain = World::getInstance().getDomain();
   for (size_t i=0;i<2;++i) {
-    for (std::vector<size_t>::const_iterator iter = state->bondside_sets[i].begin();
+    for (BoostGraphHelpers::Nodeset_t::const_iterator iter = state->bondside_sets[i].begin();
         iter != state->bondside_sets[i].end(); ++iter) {
       atom &walker = *World::getInstance().getAtom(AtomById(*iter));

src/Actions/MoleculeAction/StretchBondAction.def

@@ -21 +21 @@
 (RangeValidator< double >(0., 10.))
 
-#define statetypes (std::vector<Vector>)(std::vector< std::vector<size_t> >)(const molecule *)
+#define statetypes (std::vector<Vector>)(std::vector< BoostGraphHelpers::Nodeset_t >)(const molecule *)
 #define statereferences (Shift)(bondside_sets)(mol)