Index: src/Graph/BondGraph.cpp
===================================================================
--- src/Graph/BondGraph.cpp	(revision 326000af78937d00c935278d1a92caff13b936cd)
+++ src/Graph/BondGraph.cpp	(revision 4f04ab893d4a13da0d15780297e292e432c9536a)
@@ -33,6 +33,15 @@
 #endif
 
+// boost::graph uses placement new
+#include <boost/graph/adjacency_list.hpp>
+
 #include "CodePatterns/MemDebug.hpp"
 
+#include <algorithm>
+#include <boost/bimap.hpp>
+#include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+#include <boost/function.hpp>
+#include <boost/graph/max_cardinality_matching.hpp>
 #include <deque>
 #include <iostream>
@@ -320,2 +329,234 @@
   return treeedges;
 }
+
+struct hasDelta {
+  bool operator()(atom *_atom) {
+    const double delta =
+        _atom->getType()->getNoValenceOrbitals() - _atom->CountBonds();
+    return (fabs(delta) > 0);
+  }
+};
+
+typedef std::set<atom *> deltaatoms_t;
+typedef std::set<bond::ptr> deltaedges_t;
+
+static int gatherAllDeltaAtoms(
+    const deltaatoms_t &allatoms,
+    deltaatoms_t &deltaatoms)
+{
+  static hasDelta delta;
+  deltaatoms.clear();
+  for (deltaatoms_t::const_iterator iter = allatoms.begin();
+      iter != allatoms.end(); ++iter) {
+    if (delta(*iter))
+      deltaatoms.insert(*iter);
+  }
+  return deltaatoms.size();
+}
+
+typedef boost::bimap<int,atom*> AtomIndexLookup_t;
+
+static AtomIndexLookup_t createAtomIndexLookup(
+    const deltaedges_t &edges)
+{
+  AtomIndexLookup_t AtomIndexLookup;
+  size_t index = 0;
+  for (deltaedges_t::const_iterator edgeiter = edges.begin();
+      edgeiter != edges.end(); ++edgeiter) {
+    const bond::ptr & _bond = *edgeiter;
+
+    // insert left into lookup
+    std::pair< AtomIndexLookup_t::iterator, bool > inserter =
+        AtomIndexLookup.insert( AtomIndexLookup_t::value_type(index, _bond->leftatom) );
+    if (inserter.second)
+      ++index;
+
+    // insert right into lookup
+    inserter = AtomIndexLookup.insert( AtomIndexLookup_t::value_type(index, _bond->rightatom) );
+    if (inserter.second)
+      ++index;
+  }
+  return AtomIndexLookup;
+}
+
+typedef std::map< std::pair<atom *, atom*>, bond::ptr> BondLookup_t;
+static BondLookup_t createBondLookup(
+    const deltaedges_t &edges)
+{
+  BondLookup_t BondLookup;
+  for (deltaedges_t::const_iterator edgeiter = edges.begin();
+      edgeiter != edges.end(); ++edgeiter) {
+    const bond::ptr & _bond = *edgeiter;
+
+    // insert bond into pair lookup
+    BondLookup.insert(
+        std::make_pair(
+            std::make_pair(_bond->leftatom, _bond->rightatom),
+            _bond)
+    );
+  }
+  return BondLookup;
+}
+
+typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> graph_t;
+typedef std::vector<boost::graph_traits<graph_t>::vertex_descriptor> mate_t;
+
+static void fillEdgesIntoMatching(
+    graph_t &g,
+    mate_t &mate,
+    const AtomIndexLookup_t &AtomIndexLookup,
+    const BondLookup_t &BondLookup,
+    deltaedges_t &matching
+    )
+{
+  matching.clear();
+  boost::graph_traits<graph_t>::vertex_iterator vi, vi_end;
+  for(tie(vi,vi_end) = boost::vertices(g); vi != vi_end; ++vi)
+    if (mate[*vi] != boost::graph_traits<graph_t>::null_vertex() && *vi < mate[*vi]) {
+      atom * leftatom = AtomIndexLookup.left.at(*vi);
+      atom * rightatom = AtomIndexLookup.left.at(mate[*vi]);
+      std::pair<atom *,atom *> AtomPair(leftatom, rightatom);
+      std::pair<atom *,atom *> reverseAtomPair(rightatom, leftatom);
+      BondLookup_t::const_iterator iter = BondLookup.find(AtomPair);
+      if (iter != BondLookup.end()) {
+        matching.insert(iter->second);
+      } else {
+        iter = BondLookup.find(reverseAtomPair);
+        if (iter != BondLookup.end()) {
+          matching.insert(iter->second);
+        } else
+          ASSERT(0, "fillEdgesIntoMatching() - cannot find ("+toString(*vi)+","
+              +toString(mate[*vi])+") in BondLookup.");
+      }
+    }
+}
+
+static bool createMatching(
+    deltaedges_t &deltaedges,
+    deltaedges_t &matching
+    )
+{
+  // gather all vertices for graph in a lookup structure
+  // to translate the graphs indices to atoms and also to associated bonds
+  AtomIndexLookup_t AtomIndexLookup = createAtomIndexLookup(deltaedges);
+  BondLookup_t BondLookup = createBondLookup(deltaedges);
+  const int n_vertices = AtomIndexLookup.size();
+
+  // construct graph
+  graph_t g(n_vertices);
+
+  // add edges to graph
+  for (deltaedges_t::const_iterator edgeiter = deltaedges.begin();
+      edgeiter != deltaedges.end(); ++edgeiter) {
+    const bond::ptr & _bond = *edgeiter;
+    boost::add_edge(
+        AtomIndexLookup.right.at(_bond->leftatom),
+        AtomIndexLookup.right.at(_bond->rightatom),
+        g);
+  }
+
+  // mate structure contains unique edge partner to every vertex in matching
+  mate_t mate(n_vertices);
+
+  // get maximum matching over given edges
+  bool success = boost::checked_edmonds_maximum_cardinality_matching(g, &mate[0]);
+
+  if (success) {
+    LOG(1, "STATUS: Found a matching of size " << matching_size(g, &mate[0]) << ".");
+    fillEdgesIntoMatching(g, mate, AtomIndexLookup, BondLookup, matching);
+  } else {
+    ELOG(2, "Failed to find maximum matching.");
+  }
+
+  return success;
+}
+
+
+int BondGraph::calculateBondDegree(const deltaatoms_t &allatoms) const
+{
+  deltaatoms_t deltaatoms;
+  deltaedges_t deltaedges;
+  deltaedges_t matching;
+
+  LOG(1, "STATUS: Calculating bond degrees.");
+  if (DoLog(2)) {
+    std::stringstream output;
+    output << "INFO: All atoms are: ";
+    BOOST_FOREACH( atom *_atom, allatoms) {
+      output << *_atom << " ";
+    }
+    LOG(2, output.str());
+  }
+
+  size_t IterCount = 0;
+  // 1. gather all atoms with delta > 0
+  while ((gatherAllDeltaAtoms(allatoms, deltaatoms) != 0) && (IterCount < 100)) {
+    if (DoLog(3)) {
+      std::stringstream output;
+      output << "DEBUG: Current delta atoms are: ";
+      BOOST_FOREACH( atom *_atom, deltaatoms) {
+        output << *_atom << " ";
+      }
+      LOG(3, output.str());
+    }
+
+    // 2. gather all edges that connect atoms with both(!) having delta > 0
+    deltaedges.clear();
+    for (deltaatoms_t::const_iterator atomiter = deltaatoms.begin();
+        atomiter != deltaatoms.end(); ++atomiter) {
+      const atom * const _atom = *atomiter;
+      const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
+      for (BondList::const_iterator bonditer = ListOfBonds.begin();
+          bonditer != ListOfBonds.end();++bonditer) {
+        const bond::ptr _bond = *bonditer;
+        if ((_bond->leftatom == _atom) && (deltaatoms.find(_bond->rightatom) != deltaatoms.end()))
+          deltaedges.insert(_bond);
+        else if ((_bond->rightatom == _atom) && (deltaatoms.find(_bond->leftatom) != deltaatoms.end()))
+          deltaedges.insert(_bond);
+      }
+    }
+    if (DoLog(3)) {
+      std::stringstream output;
+      output << "DEBUG: Current delta bonds are: ";
+      BOOST_FOREACH( bond::ptr _bond, deltaedges) {
+        output << *_bond << " ";
+      }
+      LOG(3, output.str());
+    }
+    if (deltaedges.empty())
+      break;
+
+    // 3. build matching over these edges
+    if (!createMatching(deltaedges, matching))
+      break;
+    if (DoLog(2)) {
+      std::stringstream output;
+      output << "DEBUG: Resulting matching is: ";
+      BOOST_FOREACH( bond::ptr _bond, matching) {
+        output << *_bond << " ";
+      }
+      LOG(2, output.str());
+    }
+
+    // 4. increase degree for each edge of the matching
+    LOG(2, "DEBUG: Increasing bond degrees by one.");
+    for (deltaedges_t::iterator edgeiter = matching.begin();
+        edgeiter != matching.end(); ++edgeiter) {
+      (*edgeiter)->setDegree( (*edgeiter)->getDegree()+1 );
+    }
+
+    // 6. stop after a given number of iterations or when all atoms are happy.
+    ++IterCount;
+  };
+
+  // check whether we still have deltaatoms
+  {
+    hasDelta delta;
+    for (deltaatoms_t::const_iterator iter = allatoms.begin();
+        iter != allatoms.end(); ++iter)
+      if (delta(*iter))
+        ELOG(2, "Could not find satisfy charge neutrality for atom " << **iter << ".");
+  }
+
+  return deltaedges.size();
+}
Index: src/Graph/BondGraph.hpp
===================================================================
--- src/Graph/BondGraph.hpp	(revision 326000af78937d00c935278d1a92caff13b936cd)
+++ src/Graph/BondGraph.hpp	(revision 4f04ab893d4a13da0d15780297e292e432c9536a)
@@ -493,22 +493,11 @@
       AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
   {
-    const std::set<bond::ptr> backedges = getBackEdges();
-    const std::set<bond::ptr> treeedges = getTreeEdges();
-
-    //LOG(1, "Correcting Bond degree of each bond ... ");
-    int No = 0, OldNo = -1;
-    do {
-      OldNo = No;
-      No=0;
-      for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
-        No+=CorrectBondDegree(*AtomRunner, backedges);
-      }
-      for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
-        No+=CorrectBondDegree(*AtomRunner, treeedges);
-      }
-    } while (OldNo != No);
-    //LOG(0, " done.");
-    return No;
-  }
+    std::set<atom *> allatoms;
+    for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner)
+      allatoms.insert(*AtomRunner);
+    return calculateBondDegree(allatoms);
+  }
+
+  int calculateBondDegree(const std::set<atom *> &allatoms) const;
 
   bool operator==(const periodentafel &other) const;
