Index: src/Actions/MoleculeAction/StretchBondAction.cpp
===================================================================
--- src/Actions/MoleculeAction/StretchBondAction.cpp	(revision 6afd461a51752a0e6a52f5d0ae8ee1b9449deec8)
+++ src/Actions/MoleculeAction/StretchBondAction.cpp	(revision 89235eadd38fcd3984c6925d02b42e53098bbae0)
@@ -33,7 +33,4 @@
 #endif
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/breadth_first_search.hpp>
-
 //#include "CodePatterns/MemDebug.hpp"
 
@@ -52,4 +49,6 @@
 #include "Descriptors/AtomIdDescriptor.hpp"
 #include "Graph/BoostGraphCreator.hpp"
+#include "Graph/BoostGraphHelpers.hpp"
+#include "Graph/BreadthFirstSearchGatherer.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
@@ -61,32 +60,4 @@
 #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,4 +81,17 @@
     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,14 +99,11 @@
   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,64 +114,47 @@
   // 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,5 +178,5 @@
   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,5 +194,5 @@
   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));
Index: src/Actions/MoleculeAction/StretchBondAction.def
===================================================================
--- src/Actions/MoleculeAction/StretchBondAction.def	(revision 6afd461a51752a0e6a52f5d0ae8ee1b9449deec8)
+++ src/Actions/MoleculeAction/StretchBondAction.def	(revision 89235eadd38fcd3984c6925d02b42e53098bbae0)
@@ -21,5 +21,5 @@
 (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)
 
